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The Solution of Aeroelastic Problems by 
Means of Electrical Analogies 


R. H. MacNEAL,* G. D. McCANN? ann C. H. WILTS$ 
California Institute of Technology 


SUMMARY 


This paper describes a general method of mathematical anal 
The _ tech- 
nique permits the representation of both the air-frame elastic 
structure and the aerodynamic forces by electrical analogies to 
Portions of 


ysis recently developed for aeroelastic problems. 


preserve the original physical form of the problem. 
the air frame which need to be represented as distributed elastic 
structures can be set up as finite difference equations for beams in 
combined bending and torsion, or, if broad short wings are in 
volved, suitable analogies have been developed for two-dimen 
sional plates and box structures. The physical stations for each 
part of the air frame are thus preserved, and the aerodynamic 
forces or torques are introduced at the desired stations by ap- 
propriate electric circuit techniques 

In this manner the major portions of an air frame can be repre 
sented to a high degree of complexity, together with all important 
appendages such as landing gear, nacelles, and control surfaces 
One important advantage of the technique is that it not only 
permits the analysis of complex systems but also general studies of 
such systems for optimum design purposes 


(1) INTRODUCTION 


— GENERAL FIELD OF AEROELASTICITY embraces 
aircraft problems in which forces of three types 
play an essential role, these forces being aerodynamic 
forces, elastic forces, and inertia forces. In recent 
years, the number of aircraft problems in which all three 
types of forces must be considered has steadily in- 
creased. At the present time a list of such problems will 
include! (1) flutter, (2) buffeting, (3) stability and 
control, (4) divergence, (5) reversal of control, and (6) 
gusts. 

In addition to these problems, there are two other re- 
lated types, mechanical vibrations and landing im- 
pacts, in which the aerodynamic forces do not as yet 


play an important part. 
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Another recent trend points toward the interde- 
pendence of the above types of problems with the at- 
tendent need of a unified approach. To quote an 


eminent British authority, Collar:! 


‘In short, it is evident that we are no longer dealing with 
a series of subjects each in its own watertight compartment: 
there is a definite coalescence of the subjects into an inte- 
grated whole, which may be defined as the dynamics of a de 
formable airplane. And we are faced with the question: 
what are to be our methods of treatment of this unified prob 


lem?” 


As a consequence of this trend, the number of degrees 
of freedom that must be taken into account in the solu- 
tion of any one type of problem has become extremely 
large so that the use of mechanical computational aids 
is almost a sine qua non. It is the purpose of this paper 
to present a method of analog computation which it is 
believed is well suited to the solution of aeroelastic 
problems and which closely approaches the ideal of a 
unified approach to the ‘dynamics of a deformable air- 
plane.”’ 

The analogy method to be presented here is similar 
to that employed with the electronic differential analyzer 
type of computer (such as the REAC)? in that the co- 
ordinates of the problem are represented by voltages 
in the computer. It differs from the method employed 
with these computers in that advantage is taken of the 
similarity in form of the equations of a general linear 
mechanical system with a finite number of degrees of 
freedom (Lagrange’s equations) and the equations of a 
Two quan- 
namely, 


linear electrical network (Kirchoff’s laws). 
tities with fundamentally different properties 
coordinates and generalized forces—occur in a mechani- 
cal system. The voltages and currents of an electrical 
circuit have exactly the same functional relationship as 
the coordinates and forces of a mechanical system. 
Similarly, ‘impedances’ that express invariant prop 
erties of a mechanical system (masses, springs, etc.) 
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have an exact counterpart in the ‘“‘impedances’’ of an 
electrical circuit (condensers, inductors, etc.). If the 
mechanical system consists of a few simple elements, 
then it is likely that the analogous electrical circuit is 
equally simple. 

The method of computation to be described consists 
of setting up an electrical network having a direct cor- 
respondence of all its properties with those of the sys- 
tem under investigation (current-force, voltage-velocity, 
mass-capacitance, etc.) and then of measuring the cur- 
rents and voltages induced by appropriate “forcing 
functions.”’ 

In the differential analyzer type of computer, on the 
other hand, quantities of only one type (voltages) are 
recognized so that forces and coordinates must be rep- 
resented in the same way. The solution to a problem 
is obtained by electrically integrating the equations of 
motion. This integration is carried out by amplifiers 
whose output voltage is equal to the integral with re- 
spect to time of their input voltage. 

Just as the principles of operation of the two types of 
analog computers vary, so does the type of equipment 
that they employ. The direct analog type of com- 
puter consists largely of linear passive circuit elements 
(resistors, inductors, capacitors, and transformers), 
since most of the elements of a physical system can 
usually be represented by such elements. Amplifiers 
and nonlinear devices are used to represent system 
properties that are “active’’ or nonlinear in nature. 
Where a choice exists between the use of amplifiers and 
the use of passive circuit elements, the latter will usually 
be employed. 

The differential analyzer type of computer consists 
largely of amplifiers. Although passive elements are 
employed, they function as part of the amplifier cir- 
cuits. (Condensers are used in the feedback loop of 
integrating amplifiers and potentiometers are used as 
gain controls.) The application of the differential 
analyzer type of computer to the solution of aeroelastic 
problems has recently been discussed in a paper by 
Winson.* That paper and this one provide a compari- 
son of the two methods of computation. 

The authors believe that the comparative virtues of 
the direct analogy computer for the solution of aero- 
elastic problems are the following: 

(1) The physical problem can be visualized better 
because of the similarity in form between the system 
under investigation and the electrical analog. 

(2) 
be measured (as currents), as well as the coordinates. 
This is particularly true for the beam or plate analogies 
where deflection, slopes, moments, shears, and loads 


More information is provided because forces can 


can all be measured. 

(3) The elastic and inertia forces of an air frame can 
be represented by an electrical network consisting of 
condensers, inductors, and transformers only. Usually, 
the number of condensers and inductors required is no 


greater than the number of independent elastic and 
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inertia parameters defined for the system. Amplifiers 
are primarily restricted to the generation of aerody- 
namic forces. , 

The use of passive circuit elements rather than ampli- 
fiers is considered advantageous because of their smaller 
cost, greater reliability, and inherent simplicity. 

(4) Because of the inherent simplicity that results 
from the use of passive circuit elements, systems with a 
large number of degrees of freedom can be investigated, 
Elastic systems of aircraft with more than 70 degrees 
of freedom in the elastic structure alone have been 
studied on the Cal Tech Electric Analog Computer. 

(5) The time scale for the direct analog type of 
computer is usually much faster than for the other type, 
which is important when a large amount of data must 
be obtained. 

(6) Modification of the system indicated by prelim- 
inary analysis can be more rapidly accomplished. This 
rather important point is discussed in more detail below. 

An important consideration for both types of com- 
putation is the kind of coordinates used in the aero- 
elastic system under investigation. These coordinates 
may be the translational displacements and rotations 
of points within the air frame, or they may be normal 
coordinates corresponding to the normal modes of free 
The coordinates that need to be considered 
Hence, there is 


vibration. 
are ysually fewer for the second type. 
an important advantage in the use of normal modes for 
a computational method that is embarrassed by the 
presence of a large number of degrees of freedom. 
There are, however, disadvantages in the use of nor- 
mal modes. These modes are first computed, and then 
the normalized masses and spring constants appropriate 
to each mode are obtained by integration over the entire 
system. If, in the course of the analysis of the prob- 
lem, it is desired to observe the effects of variations in 
the properties of one of the elements of the airplane 
(such as the variation of nacelle location or the vari- 
ation of tip-tank weight), it becomes necessary to re- 
compute all normalized mass and spring constants for 
each case. 

This difficulty is not present when displacements and 
rotations of points within the air frame are used as co- 
ordinates because each element of the airplane in- 
fluences only a small number of quantities in the elec- 
trical network. (In the case of the variation of a tip- 
tank weight and inertia, three condensers at most are 
affected.) This type of coordinate is well suited for 
use with the direct analog type of computer. Dis- 
tributed elastic components, such as a wing with coupled 
bending and torsion, are treated by means of finite dif- 
ference equations in which the deflections and rotations 
of a finite number of points are preserved as coordinates. 

The use of the direct component type of analog com- 
puter makes possible new techniques in aeroelastic de- 
sign. If a design is unsatisfactory (let us say that the 
flutter speed is too low), then remedial measures may be 
tried on the electrical model and their effects observed. 
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THE SOLUTION OF 
If the changes made are of a local nature, the setup time 
to change from one case to the next will be short. In 
fact, optimum design parameters may be obtained in 
this way. To take a flutter analysis as a specific prob- 
lem, the combination of variable structural parameters 
may be obtained which will have a minimum overall 
weight and will still allow the airplane to meet the re- 
quired specification as to flutter speed. 

Besides this, general studies of important elastic 
problems can be made on an analog computer. It is not 
uncommon for such a study to consist of several hundred 
cases, each case involving a relatively minor modifica- 
tion from the preceding case. 

The principle disadvantages of the direct analog 
type of computer are : 

(1) It is inherently less accurate because of the 
presence of parasitic-resistance in the inductors and 
transformers. However, its accuracy is considered to 
be more than adequate for the type of problems con- 
sidered here. 

(2) The design of nonlinear devices is made more 
difficult by the faster time scale employed. This ob- 
jection is not relevant to linearized aeroelastic prob- 
lems. 

In the remainder of this paper, the application of 
direct analog computation to aeroelastic problems will 
be examined in detail, and some illustrative examples 
will be considered. 

ELASTIC 


(II) EvLectrricAL ANALOGIES FOR’ THE 


FORCES AND INERTIA FORCES 


The similarity in form between the equations of a 
mechanical system and the equations of an electrical 
network was mentioned in the introduction to this 
paper. Two distinct electric circuits may be derived 
depending on whether the Lagrangian equations of the 
mechanical system are made equivalent to Kirchoff's 
voltage law (the sum of voltages in each loop = OQ) or 
Kirchoff's current law (the sum of currents entering 
each node = Q). 

The analogous relationships that result from these 
two methods are summarized in Table 1. 

The second type of analogy in which force is analogous 
to current, the so-called ‘‘mobility analogy,”’ has been 
used exclusively by the authors in aeroelastic applica- 
tions. When damping is neglected in an elastic system 


TABLE | 
Electrical-Mechanical Analogies 
Analogous Electrical Quantity 


From Kirchoff's From Kirchoff's 


Mechanical 
current law 


Quantity voltage law 
Velocity Current Voltage 
Force Voltage Current 
Mass Inductance Capacitance 


Reciprocal spring stiff- 
ness (or compliance ) 
Damping coefficient 


Inductance 
Reciprocal Resistance 
(or conductance ) 

Transformer 


Capacitance 
Resistance 


Lever Transformer 
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Fo(t) 


Loe 
ae ve 
Ke 
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M, 
& 
Ky 
MECHANICAL SYSTEM 
Fic. 1. 


Te(t) 











m 


ELECTRICAL ANALOGY 


Ve Keo M,* C 
V2=e2 keris Ma* Co 
Fo*Ig 


Simple analogous circuit. 


(which can usually be done for air frames), its electrical 
analogy contains no resistors. 

The method of setting up an electrical analogy will 
be demonstrated with reference to Fig. 1. In that 
figure the equations of equilibrium for the mechanical 


system are 


Mo + (A, + x) f V1 dt — Ke / Ve dt = ( (1) 
at . 


-K: f Vv} dt + M _ + K. f Ve dt = F,(t) (2 
€ 


The coordinates v; and v2 are velocities. In addition 
to these equations, appropriate initial conditions must 
be specified. The equations for the sum of the currents 


entering each node of the electrical circuit are‘ 


as (; yf 3 [ 4 
( e, dt — edt=Q0 (3) 
ra ee | eee 


| des l 
i e, dt + C. — fea- I,(t) (4) 
A ws “u* L 


[.(t) is a current produced by a generator with an in- 
finite source impedance. By comparing the two sets of 


equations, we see that they are identical, provided that 


v7 = @ Ky = | Ly M, = Ci 
Ke = l Le Mp = C, 


Yo = € : 


Fy, = Is 
If these conditions are satisfied, the two systems are 
analogous. In the Appendix, these conditions are 
modified to permit the introduction of scale factors. 
In accordance with established notation, d/dt is re- 
placed by the symbol p and JS dt is replaced by the 
symbol 1 /p in Eqs. (3) and (4), but the introduction of 
these symbols does not imply that the variables are re- 


placed by their Laplace transforms. 


l l l 
, — _ = 5 
(« ip + Lip + ip Lop ™ c 
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MECHANIC AL SYSTEM 
quanity | pire. FINITE DIFFERENCE EQUATIONS ee MSRY 
EQNS DISPLACEMENTS VELOCITIES CIRCUIT EQUATIONS 
AS COORDINATES AS COORDINATES 
DEFLECTION h h ph En 
z Oh 1, -h — (ph)ngi (Phin —_ | 
SLOPE (1) =e Gorge “P| (06 n+g + (PRIn) {Phin ® (Eanes © Eninei(En)n (Eno 
38 oe LU ‘) | 
MOMENT (2) | M-EI5S | Mn=-Ki(@nty-On-g) | Mas “EL[ipone $ -(P@m-5] | lighy = — [SPInee = Eades | 
SHEAR (3) Qs 37 Qn+5= Mos Mn Qn+5* Mosi= Mp (i, ne¢* Cones (pn 
Q “is if ; ; ? : 
LOAD (4) ~~ Qn v~ Cnth Gor Ayan =- [Qn+4-On-¢ ] liq = — [inn +g - linn] 
INERTIA s-mazh =-m 2h --m A(ph er 
LOAD (5) iil hated qe—m 5. ) (ign =-c STA 
ANALOGOUS QUANTITIES: 
A 
En=ph Zslp ETP 
Eg = pe C = mAy 
in = Q T * Ay (TRANSFORMER TURNS RATIO) 
ie =M p = 4 
ig = Ay.q 


Table 2. 


I 
= Lap € + (ce + p)e = Jo(f) (6) 


The mechanical system in Fig. | can also be visualized 
as a torsional system in which v, and v2 represent angu- 
lar velocities. This torsional system can be thought of 
as a representation by two finite-difference cells of a dis- 
tributed shaft vibrating in torsion, one end of which is 


clamped. In this manner it can be shown that the elec- 
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Fic. 2. Electrical analogy for the bending of a beam. 


Development of the finite difference analogy for the bending of a beam 


trical analogy for a shaft vibrating in torsion is a lumped 
transmission line with series inductors and shunt ca 
pacitors. 

The next example to be considered is the finite dif 
ference analogy for the bending deflections of a beam. 
This analogy is described in more detail in reference 
». Itis of utmost importance for the direct analog solu- 
tion of aeroelastic problems, since the major structural 
components of an airplane (wing, fuselage, and tail as- 
sembly) all have flexural stiffness. The derivation of 
this analogy will be demonstrated with reference to 
Fig. 2 and Table 2. 

In Table 2, the transformation from differential equa 
tions to finite difference equations in which velocities 
are used as coordinates is shown in the first three col- 
umns. A mechanical system that satisfies the finite 
difference equations is shown in Fig. 2. Weightless 
levers are used to transform the difference of the deflec- 
tions of adjacent pivot points into the deflections of 
a set of springs that represent the bending stiffness of 
the beam. Dynamic loads are represented by a set 
of lumped masses. 

The equations of the electrical analogy are obtained 
by the substitutions listed in Table 2. It can easily be 
shown that the electrical circuit satisfies these equa- 
tions. The ratio of the voltages across the two coils of 
an ideal transformer is proportional to the ratio of the 
number of turns in the two coils. If the turns ratio 7 
is taken equal to Ay as shown, the voltages on the two 
sides of the transformer satisfy the first equation. The 
second equation is Faraday’s Law of Induction for the 
The third 


current in the inductors of the /, circuit. 
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equation relates the currents in the two coils of the 
Since currents transform inversely 
The 
fourth equation is Kirchoff’s current law for the nth 


ideal transformer. 
as the turns ratio, this equation is also satisfied. 


node, while the last equation defines the relationship 
between current and voltage of a capacitor. From this 
derivation it is seen that the electrical analogy for the 
lever is the transformer. 

The method of applying boundary conditions to the 
beam analogy is illustrated in Fig. 3, which shows a 
five-cell representation of a cantilever beam with uni 
form mass loading. The boundary conditions at the 
clamped end of the beam are that the deflection and 
slope are both equal to zero; since these quantities are 
respectively analogous to the time integrals of the volt 
ages :, and /, the boundary conditions are satisfied 
by shorting these voltages to ground at the nodes cor- 
responding to the location of the clamped end. The 
boundary conditions at the free end are that the bend 
ing moment and shear are both equal to zero; since these 
quantities are respectively analogous to the currents 
i, and 7,, these boundary conditions are satisfied by 
opening the slope and deflection circuits in the loops 
corresponding to the location of the free end. 

There are many interesting topics that have arisen 
in connection with beam analogies for which the avail- 
able space does not permit an adequate discussion. 
Some of these topics are: 

(1) The development of new analogies that have a 
smaller error due to finite difference approximation. 
One such analogy, which represents the elastic forces in 
an exact manner, has been derived by Russell. 

(2) The representation of additional phenomena 
such as the effects of shear,’ rotary inertia,° axial stress, 
beam curvature,® and unequal lumping.® 

(3) The investigation of parasitic effects such as 
transformer leakage and magnetizing inductance® and 
the series resistance of transformers and inductors. 

(4) The extension of the beam analogy to the theory 
of thin plates’ and the further extension to thick plates 
and shells. 

The application of the beam analogy to an aircraft 
vibration problem is illustrated in Fig. 4, which repre- 
sents a wing vibrating in bending and torsion with in- 
In this case, the fuselage is represented 
Mass coupling due to the 


ertia coupling. 
by its inertia forces only. 
first moment S is represented by the condensers C,,,. 
The first moment is negative for section 1. This is ac- 
complished through use of a transformer that reverses 
the sign of the voltage E, = pa. This illustrates a 
second important use of the transformer, which is to 
eliminate negative circuit elements that would other- 
wise arise in connection with coupling terms. When 
scale factors are introduced, the circuit elements will 
be changed in the manner given in the Appendix. 

A third important application of the transformer is 
its use in the transformation of coordinates, particularly 
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Analogy for the vibration of an airplane wing 
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Transformation from elastic axis coordinates to 
air-stream coordinates. 


Fic. 5. 


in the rotation of coordinate axes. This will be illus 
trated with reference to Fig. 5, which shows an aircraft 
with a swept wing. If the wing is treated as a beam, it 
is most convenient to use elastic axis coordinates. The 
axis of twist a,, is the elastic axis, and the axis of slope 
rotation @ is perpendicular to the elastic axes. In order 
to attach the wing to the fuselage or to introduce aero- 
dynamic forces, these coordinates must be transformed 
into the pitch angle a and the roll angle 0. If u is the 
angle of sweepback, the equations expressing this trans 


formation are 


a =a,,cosu+ @6sin pu ( 
0 = 


St 


—a,, sin wp + 6 cos u ( 


¢a 
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Fic. 7. Electrical circuits for the aerodynamic lag function 


For more convenient electrical representation, these 


equations are divided by cos u. 
= a, + 0 tan pu (9) 
—a,, tan p + 0 (10) 


a/COS b 


3/cos wp = 


The transformer network shown in Fig. 5 carries out 
this transformation by the simple addition of voltages. 
The bending and twisting moments are also correctly 
transformed which can be proved for a general coordi- 
nate transformation by using the knowledge that no 
‘energy is stored or dissipated in an ideal transformer.* 

(III) ExvecrricaL ANALOGIES FOR THE AERODYNAMIC 

FORCES 


In the usual small motion aerodynamic theory, the 
aerodynamic forces are linear in the displacements of 
the airfoil. These forces, however, are not, in general, 
conservative, passive, or “bilateral’’ in the sense that 
the matrix of force coefficients is symmetrical with re- 
spect to the principal diagonal. Hence, the electrical 
analogy for these forces must include resistors and 
amplifiers. 

In the discussion that follows, those terms that are 
proportional to acceleration of the coordinates (the 
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aero-inertia terms) have been deleted, since they can be 
represented by a system of condensers which can be 
combined with the condensers simulating the inertia of 
the air frame. 
the lifting surface is divided into strips parallel to the 


The usual assumption is made that, if 


air stream, then the forces on each strip can be com- 
puted independently. The possibility of relaxing this 
assumption is briefly discussed at the end of Part ITI, 
The coordinates that are used are measured with refer- 
ence to the initial undisturbed position of the airfoil. 
Subject to these qualifications, the equations for the un- 
steady lift and moment on the strip shown in Fig. 6, 


moving in an incompressible fluid, are*~'! 


. i cp\ {U c 
~e (SP) Tjyv? * 
.f oP ph X2 )! 
(P(r + pa + 0 (11 
M.. = —mL —- 2ra( “ie 
SV. 


The symbols in these equations which are not shown in 


L= 


Fig. 6 have the following meanings: 


q = (1/2) pV? 
U’ = gust velocity 
V = air-stream velocity 


= lag function for gusts 


= Theodorsen's function 


5 = area of strip 


Strictly speaking, the above equations are equations 
for the Laplace transforms!‘ of the variables. The use 
of operational methods in the airfoil theory has been 
The function 

These func- 


demonstrated in a paper by Sears.'” 
C(cp/2V) is there designated 1 — #(p). 
tions can be written in the terms of the modified Bessell 
function Ay and A,, but, in order that they may be rep- 
resented by electrical circuits, these functions must be 
approximated by the ratio of two polynomials in /p. 
The quantity (ph/V) + (%2/V)pa + a is the angle of 
attack at the three-quarter chord. 

Before introducing the electrical circuits, Eqs. (11) 
and (12) must be transformed into the equations for the 


analogous electrical quantities. In accordance with 
the theory developed in the Appendix, we take 

b =p/N pa = (—Kb/N)E, 

ph = (Ka/N)E, M.. (—K/b)i, 

L = (K/a)t, 
where V is a time-scale factor (mechanical time = .V X 


electrical time) and A is an arbitrary dimensionless 
scale factor relating the energy levels of the two systems. 
The air-stream velocity |’ may be eliminated as a vari- 


able parameter by choosing 
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a NoVo _ alo by Vo (13) and (14) are represented by a functional block dia- 
a a ilies oe co gram in Fig. Sa. This diagram may be thought of as a 
tcpological representation of the operations indicated 
where |’) is the “reference” velocity. By meansof these jy, Eqs. (13) and (14) rather than as an electrical cir- 


substitutions, Eqs. (11) and (12) become cuit. The block labeled —x, — (.No\o/p) has this ex- 
. aobo § cp No 7a pression as its transfer function. The blocks labeled 
in = 24 qQod VoVo ) -u(; <P _)(*) U + 4 i + G, and G, produce the currents 7, and /). 

An electrical circuit performing these operations in an 


‘( cP \(- Ho =o” es VoVo E )! (13) identical manner is shown in Fig. Sb. The transfer 
bo - a a 


IN y - . y , . ° ° 
2NoVo Pp | function —x. — (NoVo/p) is accomplished by using an 
boxy . — be? \ fe? amplifier G, to generate a current proportional to the 
—™ de i, — <8ge NoVo/\8 Fa (14) voltage ,. The voltage across the series combination 
of elements R,; and C, is this current times the imped- 
where go = (1/2)pVo?. The lift and moment currents — ance of combination 


are now independent of velocity, but the inductors rep- 
resenting elastic forces have become functions of ve-  — 








locity, since, by the Appendix, 























tr 

i, ik , V0" - —_ [ | 

..* a;'a;'Kiy = (a;')o(aj’)o(K iy) [7 (15) ine P cup) 

“ij | ee 4 

The condensers representing inertia forces are still 
independent of velocity. Eh vee) 
] | 95) | 
Cis = [(a;')o(aj’)o/ No?) Wis (16) tn 

i. bof, c® 6 ° , Pe 
Re" 2*%S Ry,” & 6, * 2va,S. ee Ge" BE+X,Gp 


The above discussion illustrates the interesting 


° : : : 6c FUNCTIONAL DIAGRA Ri R 
theorem that an increase in the velocity of an elastic BLOCK DIAGRAM OF AERODYNAMIC CIRCUIT 








body, moving in an incompressible fluid, is equivalent 
to a decrease in the stiffness. Electrically, it is easier to 
allow the elastic inductors to vary with speed than it is 
to alter the aerodynamic circuits. 

Theodorsen’s function C may be accurately repre- 





sented by an expression of the following type: 


(1 t (] ts 
gees (17) 
(1 + top) (1 + tap) 








as i RR; a a: ee See eee ra 8b ELECTRICAL CIRCUIT THAT IS FUNCTIONALLY IDENTICAL WITH THE 
In a paper by Biot and Wianko,"* a simple electrical net eaten aI a 


work that has this transfer function is derived and used 
in the solution of the flutter problem. Their work is 
probably the earliest using the direct electrical analogy 
method in such problems. Theirs or more complex 
forms can be represented by linear transfer functions 
and, hence, by electrical networks. A circuit for the 
transfer function of Eq. (17) is shown in Fig. 7a. For 
illustrative purposes we shall show a simpler approxi- 
mation to C(p) which is frequently sufficient for the 
solution of subsonic problems. This transfer function 
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@c_ACTUAL CIRCUIT USED 


is Fic. 8. Circuit for generating aerodynamic forces 





1+ hp 


C(p) = (f2 = 2t) (1S) TO WING STATIONS 


1 + bop Ri 


which, for imaginary values of / = iw, is a semicircle 


7 


in the complex plane with center at x = 0.75, y = 0, 1 
and radius = 0.25. The electrical circuit, together T if : ® | | | | ; 
with circuit element values, is shown in Fig. 7b. These l a 
networks are designed to be fed by a voltage source — SHagE e0ceD [¥H ae ee ee 

i , J i GUST C(B)] MEAN CHORD 
with zero output impedance and to feed into networks ane , 

‘ 5 ll é ‘ , s Fic. 9. Circuit for generating angle of attack due to a sharp- 
with infinite (or very high) input impedance. Eqs. edged gust 
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FE, = [R, + (1 Cip) | (—G,)(-,) (19) 


RC, is equal to x2/ No Vo. 
Addition is performed by means of the resistors R, 
and R,,. 


channel of such an adding circuit is inversely propor- 


It may easily be shown that the gain of any 


tional to the resistance in the channel, independently 
of the number of channels or the input impedance of the 


circuit following the adder. Hence, 


A,/Aq = R,/R, (20) 


The values of such adding resistors are always large. 
The simple approximation of Fig. 7b to the function 
C(p) is accomplished by R2, C:, and R;. The amplifier 
(—A,;) is used to achieve a change in sign of the un- 
lagged term before entering the second adder. The 
operation of the elements Ry, G,, and G, is identical 
with their operation in Fig. Sa. 

In line with the policy of a minimum use of amplifiers, 
that of Fig. Sc 
This is accomplished by 


the actual circuit used contains only 
three amplifiers per cell. 
eliminating the second adding circuit of Fig. Sb and re- 


placing the current generator G, by a transformer with 
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a turns ratio 7 = G,/G, = boxi/do. The function yp) + 
C(p) by which the gust velocity is now multiplied can 
be adequately represented by a simple circuit. The 
elements 1; and R, serve to represent the fraction of 
the unlagged /, term taken away by the high-frequency 
attenuation of the lag circuit. These elements will 
ordinarily be small. 


The circuits that have been described must be dupli- 
cated for every cell in the finite difference representa- 
tion of a wing. In the problem of gust entry for a 
straight wing of uniform chord, the function ¥(p) C(p) 
is equal for every cell, and, hence, a single forcing func- 
If the 
wing is swept back, then the delay time in the arrival 


tion may be used to generate the gust forces. 


of the gust front at each cell can be represented by an 
¥(p)/C(p 


represents a lag of relatively short duration and, for 


electrical delay line as shown in Fig. 9. 


wings of moderate taper, can probably be given values 
corresponding to the mean aerodynamic chord, making 
possible the use of the circuit in Fig. 9 for tapered 
wings. Otherwise, ¥(p) C(p) must be duplicated at 


each cell. 
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POSITIVE DIRECTIONS FOR 
BENDING AND TORSIONAL DIRECTIONS 


WING 
BEAM IN VERTICAL BENDING 
COUPLED WITH TORSION 


Fic. 10 


FUSELAGE 


BEAM IN SIDEWISE 
BENDING COUPLED WITH TORSION 


Electrical analogy for unsymmetrical modes 


CIRCUIT FOR RIGID BODY 
YAWING OF VERTICAL 
STABILIZER AND SPRING ___ 
COUPLING TO ITS TORSION 
CIRCUIT 















VERTICAL STABILIZER 


BEAM IN SIDEWISE BENDING 
COUPLED WITH TORSION 


HORIZONTAL STABILIZER 


BEAM IN VERTICAL BENDING 
COUPLED WITH TORSION 





Vibration analysis of an airplane 
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The aerodynamic theory used in deriving the elec- 
trical circuits is correct only for a wing of infinite aspect 
ratio moving in an incompressible fluid. Important 
considerations are the extent to which these conditions 
can be violated without necessitating a radical modi- 
fication of the electrical circuits and the form that the 
electrical circuits will take for low aspect ratios and 
high Mach Numbers. 

The following quotations from aerodynamicists are 
relevant to the first consideration: 


“The main conclusion to be derived from study of the nu 
merical flutter calculations is that the effect of compressibility 
on the flutter speed (wing, bending, wing torsion, no aileron ) 
for subsonic speeds with no shocks, although complicated, is 
relatively small in the usual cases and for Mach Number of 
0.7 can be allowed for by corrections of small order to the in 
compressible case results.""!! 

“The authors feel that if it is desired to obtain theoretical 
flutter speeds which are within say 10 per cent of the actual 
flutter speeds, then in most cases it will be found necessary 


to incorporate the aerodynamic span effect in the analysis.”’® 


In the latter paper, aspect ratios of 3 and 6 were con- 
sidered. The results were presented as corrections to 
C(iw) at different spanwise locations. 
appeared to be nearly the same for different types of 


These corrections 


motion leading to the inference that they could be rep- 
resented electrically by a spanwise correction to the lag 
function without interaction among neighboring cells. 


To date compressibility corrections at the Analysis 
Laboratory have been made by allowing the coefficients 
of Eqs. (11) and (12) (including the coefficients in the 
lag functions) to be functions of Mach Number as ob- 
tained from wind-tunnel tests on rigid models. It is 
possible that this procedure, together with the intro- 
duction of simple interactions (in which the lift at one 
station contains terms proportional to the angle of 
attack at adjacent stations each multiplied by a lag 
function), may be able to give useful results for most 


wing shapes at any Mach Number. 


For swept wings of low aspect ratio, the assumption 
that the elastic properties of the wing can be represented 
by those of a simple beam is no longer justified since wing 
cambering becomes important. The most direct ap- 
proach would be to treat the wing as an anisotropic 
plate,’ with aerodynamic forces represented by a com- 
plete set of interactions. At the present time such an 
approach, even if formulated, would exceed the ca- 
pacity of existing analog computers, although construc- 
tion of a computer with adequate capacity seems to be 


practically, as well as logically, possible. 


The discussion of this section has dealt only with 
wings without flaps. It is not difficult to include aero- 
dynamic forces due to ailerons, ete., and this has been 
done in the solution of several plane studies conducted 
at the Analysis Laboratory, one of which will be dis- 


cussed in Part IV. 
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PROBLEMS 


(IV) REPRESENTATIVE PROBLEMS 


DISCUSSION OF 


The solutions of three problems have been chosen for 
discussion. In the space available it is hardly possible 
to give more than a brief impression of what the prob- 
lems were and how they were solved on the Cal Tech 
Electric Analog Computer.'® 
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(a) RELATIVE ANGULAR MOTION OF CONTROL SURFACE 
PLUS TIMING WAVE. 





Fic. 18. Typical oscillograms from the elastic missile problem 
(a, b, and ec do not refer to the same case ) 


Problem I. Vibration Modes of an Airplane 


This is a problem that can be solved without the use 
of amplifiers or resistors, since it involves only elastic 
and inertia forces. The elastic deformations that were 
considered are: the vertical bending and torsion of the 
wing; the vertical and horizontal bending and torsion 
of the fuselage; the bending and torsion of the hori 
zontal stabilizer; bending and torsion of the vertical 
stabilizer and its rigid body motion in a vertical plane; 
and the vertical and horizontal motion of the engine 
elastically coupled to the fuselage. 

The complete circuit diagram as set up on the com- 
puter for unsym:netrical motion is shown in Fig. 10. 
Transformers are used in the coupling between the 
fuselage and the other members. This is a result of 
giving the scale factors a,’ (see Appendix) different 
values in the wing, fuselage, and tail; this was done so 
that the same range of values of the circuit elements 
could be used to represent the wing and stabilizer (where 
the inertia and elastic forces are small) as were used to 


represent the fuselage 
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Modes of free vibration were excited in the electrica] 
circuit by connecting a variable frequency voltage source 
(corresponding to a sinusoidal velocity of fixed ampli- 
tude) to a point in the circuit where large displacements 
were expected. The frequency of the voltage source was 
varied until the criterion for the presence of a normal 
mode was satisfied —that is, until the input current 
to the circuit (force applied to the airplane) was a mini- 
mum. The proximity of this minimum to zero current 
was a good indication of the ease with which the mode 
could be excited. At each minimum of current, the 
frequency of vibration and the voltages at all nodes 
were recorded. A total of 17 antisymmetrical modes 
were obtained in this way, and 15 others were ob 
tained from a similar circuit for the symmetrical 
vibrations. The results were later compared with vi 
bration tests on the actual airplane, and the comparison 
proved to be satisfactory. 

Single-Line Diagrams for Elastic Circuit Analogies. 
With circuits of this complexity, it becomes difficult 
to visualize the analogy in terms of the actual mechani- 
cal system. It is therefore desirable to use simplified 
single-line or symbolic circuit diagrams. Such a simpli 
fied representation is illustrated in Fig. 11 for an air 
plane in symmetric motion. Each type of coordinate 
(p60, ph, ete.) is represented by a transmission wire. 
For each elastic cell, the mass, spring, and transformer 
couplings between coordinates is symbolized as shown. 
For the junctions between individual beams, the proper 
boundary conditions are indicated by the connections 
of the coordinate wires. For instance, at the wing root 
(Fig. 11), the bending slope (or roll angle) for the wing 
is zero for symmetrical motions. Thus, the (p@) co- 
ordinate wire of the wing is grounded. The wing tor- 
sional displacement (or pitch angle) is made equal to the 
fuselage bending slope by connecting the (pa) wire of the 
wing to the (p@) wire of the fuselage. 

This method of symbolization will be used for all 


subsequent examples. 


Problem II. The Effect of Fuselage Elasticity on the 

Longitudinal Stability of a Missile 

In the past few years it has become increasingly ap- 
parent that elasticity plays an important role in the 
stability and control of all types of aircraft. The pres- 
ent study was undertaken for a Government agency 
to determine the amount by which the damping 
and frequency of the main symmetrical pitching motion 
of a missile are changed because of elastic deformation 
of the fuselage and also to determine the change in the 
effectiveness of aerodynamic surfaces to alter the flight 
path. 

A simplified schematic diagram of the electric cir- 
cuit is shown in Fig. 12. The circuits used for the aero- 
dynamic forces at the control surface, wing, and tail 
were similar to that of Fig. Sc. Several supersonic 


speeds were considered. In these cases the aerody 
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namic forces were considered to have the same form as 
for incompressible subsonic flight, except that the co 
efficients of Eqs. (11) and (12) were varied with Mach 
Number in what was considered by the aircraft engi 
neers to be an appropriate manner. The lift-slope 
coeflicient and the center of pressure were obtained 
from wind-tunnel tests, while the lag function coeffi 
cients were obtained from a paper by Stewart and Li." 
Thus, the values of the electrical circuit elements had 
to be changed in going from one speed to another. The 
linear dimensions of the aerodynamic surfaces were 
small compared to the length of the missile, and their 
rigidity was so large that it seemed reasonable to as 
sume these surfaces to have infinite rigidity. 

Downwash from control surface to wing and from wing 
to tail was accomplished by adding a term to the angle 
of attack proportional to the lift on the preceding sur 
face. The downwash delay time was simulated by an 
electrical delay line. 

The fuselage was represented by the electric analog 
of a free-free beam with six cells; a preliminary vibra 
tion mode analysis showed good agreement with vibra 
tion tests conducted on the missile. 

Transient oscillation of the missile was excited by a 
pulsed angular displacement of the control surface. 
Typical oscillograms of the transient response are 
shown in Fig. 13. Damping and frequency were ob 
tained by a numerical analysis of these records and also 
from a “data-analyzing network,’ which operated di- 
rectly on voltages from the computer. 


Problem III. Flutter Analysis of a Swept-Wing Airplane 


The new design techniques mentioned at the end of 
Part I were fully utilized in connection with this prob- 
lem in that several structural components were varied 
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Fic. 14. Block diagram of circuit for flutter model of swept-wing 
aircraft 


within the limits imposed by other design considerations 
in search of an optimum flutter design (i.e., minimum 
weight for a given flutter speed ). 

A simplified schematic diagram of the airplane is 
shown in Fig. 14. The fuselage was kept rigid but was 
free to roll, pitch, and plunge. The wing was divided 
into five cells, both for the finite difference beam an- 
alogy and for the aerodynamic loads. Since the wing 
was swept back, a transformation of the slope and twist 
of the elastic axis into air-stream coordinates as dis- 
cussed at the end of Part II was required. Additional 
features included: a nacelle with a flexible pylon to be 
attached somewhere between section | and section 2; a 
flexibly attached wing tank; and an aileron with vari 
able mass balance and control cable stiffness. The aero 
dynamic theory used was that for an incompressible 
fluid and infinite aspect ratio. Compressibility and 
aspect ratio corrections were added later. 

The aerodynamic cells for sections 0, |, and 2 were 
similar to those of Fig. Sc; but the two outboard sec- 


tions contained a number of interconnections with the 
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aileron, so that in the preliminary testing these three 
components were treated as a single unit. These tests 
consisted of calculating the flutter speed of each aero- 
dynamic unit when connected to a circuit representing 
an infinite airfoil with simple distributed spring sup- 
ports. The speed was varied by changing the induc- 
tors corresponding to these springs as explained in Part 
III. For each speed the transient decay rate of the 
free oscillation was recorded. The flutter speed was 
obtained by drawing a curve of damping vs. speed and 
by noting the point at which the curve crossed the speed 
axis. 

The same method was used for obtaining the flutter 
speed of the airplane as a whole. In this case the decay 
rate of the Jeast damped component of free oscillation 
was measured. Oscillations were excited by the ap- 
plication of a step lift to any one of the cells, the point 
of application being chosen so as to give a maximum 
amplitude to the least damped component. Both 
symmetrical and antisymmetrical cases were investi- 
pated. 

In the course of the problem, studies were made of 
the effect on flutter speed of varying the following 
things: (1) aileron mass balance, (2) nacelle pylon stiff- 
ness, (3) nacelle location, (4) wing-tank weight and 
pylon stiffness, (5) distribution of torsional material, 
and (6) aileron cable stiffness. 

Flutters of two types were encountered, these being 
classified roughly as wing flutter and aileron flutter. 
It was found that the effects of the items on the above 
list could not, in general, be considered independently ; 
in other words, the change in flutter speed for the com- 
bination of variations A plus B was not necessarily 
equal to the sum of the changes in flutter speed due to 
the variations A and B considered separately. Hence, 
the selection of the optimum group of values of the vari- 
able parameters could only be achieved after a large 
number of combinations had been investigated. Nearly 
1,000 points were obtained for the two airplanes studied 
(each point being the transient decay rate for one com- 
bination of altitude, speed, and structural configu- 
ration). Fig. 15 shows some typical oscillograms. 
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APPENDIX 


The differential equations of a linear time-dependent 
electrical circuit can be written in the following tensor 
form: 


G@e=1,...,8) G-1) 


The components of the vector e’ are a set of linearly 
independent node pair voltages. /; are currents that 
do not depend on the voltages e’ but come instead from 
“external” sources. Each component of the admit 
tance tensor };, expresses as a differential operator the 
relationship between the voltage e’? and the portion of 
the current entering the 7th node due to e’. The sys- 
tem is not necessarily conservative even if all J; = 0. 
)’,, is not necessarily equal to Y;;;. 

The equations of a linear mechanical system can be 
expressed in an exactly similar form. Eqs. (A-1) be 
come the equations of such a system if the components 
of e’/ represent a set of independent coordinate velocities, 
and /; are “external” forces. 

We wish to express the quantities involved in Eqs. 
(A-1) in terms of those of an ‘‘analogous’’ system whose 


equations can be written in the form 


> Vie - lg = ¢ 


CM al eer | (A-2) 


In order to do this we must make use of three 


analogy postulates as follows: 
(1) The independent variable (time) is linearly re 


lated in two analogous systems 


t= to + Ni (A-3) 
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THE SOLUTION OF 
(2) The total energy supplied from external sources 
is linearly related in two analogous systems 
E(t) = ky + K°E(d) (A-4) 
(3) An analogy matrix a,' can be given such that 
~k . - 
eé= Da? (i,k =1,...,2) (A-5) 
P 
As a consequence of the first two postulates, the 
power supplied from external sources is proportional in 
the two systems 
dk Idk K*dk kK 


P=—=—--= - = —p 
dj Nd Nd WN 


(A-6) 


With the aid of the above postulates, the manner in 
which the quantities /; and };; transform can be estab- 
recognized external 


lished. The power supplied by 


sources is, simply, 


P = } ie'l, (A-7) 
t 
Hence, 
> el, = (K?2/N)> eT, (A-8) 
t ? 

For our present purposes, the matrix a,' will be pre- 
sumed to contain diagonal elements a;' only. Hence, 
there is a one-to-one correspondence between the de- 
pendent variables of the two systems. 

Substitute into Eq. (A-S) for /; from Eqs. (A-1) and 
(A-2) 


dV e'Y, e? = (K? N) MeV (A-9) 
1 J 
Substitute for e' and e’ from Eqs. (A-5) 
(A-10) 


= a,;‘a; Y ,,6°6? = 


(K2/N) ¥F,,6 
ij ij 


Since the set of variables @’ is linearly dependent 


V,, = (N/K°)a;'a/'V;, (A-11) 


Also 
(A-12) 


— N ef 
> Y ,,é7 = Data vyl ) 
= (N/K) doa; ¥ ie! 
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Hence, 


I, = (N/K?®)a;'‘1, (A-13) 
Eq. (A-12) holds independently for each term of the 
Hence, all of the terms of Eqs. (A-1) 
These results can 


summation on j. 
transform in the same manner as /;. 
be summarized in the following rules: 


Given Then 
t =t +N Vij = (N/K®)aj‘ay’ Vi 
E=E,+ K°E I; = (K2/N)1/a;'; 
e = a;'é' 


A slightly different symbolism, which emphasizes 
the role of A as a dimensionless scale factor, can be ob- 
tained by replacing a,’ by (A N)a;,'.. The above rules 


become 


Given Then 
t =h+Ni Vii = (1/N)a;‘a/p Y; 
E=E,+ K°E I; = (K/a;')I 


e° = (K/N)a;‘é 
These are the transformation rules that have been 
used in the text. 
The inertia force of a generalized mass can be written 


F, = M,,(dv*/dt) 

Hence, the ‘‘admittance”’ of the generalized mass is 
Y,, = M,,;(d/dt) 

Analogously, the admittance of an electrical condenser 

is 

C;,;(N) (d/dt) 


V,, = C;,(d/di) = 


Applying the transformation rule for admittances, 


a;'a;’ 
—— 1M; 
N? 
The inductance L;; corresponding to a spring constant 
K,; can be calculated in a similar manner. 


"i 
_ ie 











Analysis of the Elastic and Plastic Stability 
of Sandwich Plates by the Method of Split 
Rigidities—IT’ 

P. P. BISLAARDt 
Cornell University 


SUMMARY 

Formulas used in the first part of this paper are derived in de- 
tail. First, the derivations are given of the reduction coefficients 
for plasticity for Cases 0 (buckling of single faces) and 1 (buckling 
of sandwich plates without taking account of shear deformation 
of core), for several cases of loading and boundary conditions. 
These reduction coefficients refer to the ratios of plastic and 
elastic buckling stresses for a given half wave length, as needed 
for the computation of sandwich plates. Using these expres- 
sions, some additional formulas are given for homogeneous plates, 
referring to the ratios of plastic and elastic buckling stresses, both 
computed for that half wave length that makes them a minimum 
Also the buckling coefficients a of Case 2 (assuming only shear 
deformation of the core to occur) are derived, as far as not yet 
‘ published elsewhere. In the Appendix formulas are given for 

application of the method to plates with anisotropic core. 


INTRODUCTION 


rr THIS SECOND PART the derivations will be given of 
those formulas used in Part I which had not yet 
been published in the author's earlier papers. A some- 
what more detailed derivation of these formulas was 
given in reference 1. The formulas involve the elastic 
buckling coefficients k) and the reduction coefficients 
for plasticity n, both referring to Cases 0 and | and ap- 
pearing in Eqs. (18) and (19) of Part I, as well as the 
buckling coefficients a for Case 2, appearing in Eq. (30) 
of Part I. To avoid misunderstandings, the equations, 
figures, and table of this paper will be numbered con- 
tinuously with those of Part I. 

References for coefficients ko, which follow from the 
theory of elastic stability, were given in Part I. A 
check of Eq. (25) for (ko), will be given near the end of 
this paper. Hence, first, the formulas for the reduc- 
tion coefficients for plasticity 7 will be derived. Then, 
Eq. (25) for (ko), will be checked. Finally, the neces- 
sary derivations of coefficients a in Eq. (30) for Case 2 
will be given. 


THE REDUCTION COEFFICIENTS FOR PLASTICITY FOR 
CASES 0 AND | 


Using the method of split rigidities, the reduction 
coefficient for plasticity refers to the ratio between the 


* Presented at the Structures Session, Eighteenth Annual Meet- 
ing, I.A.S., New York, January 23~26, 1950. The first part of 
this paper appeared in the issue of May, 1951. This second part 
contains the theoretical derivations, as far as not available else- 
where. 

+ Professor of Structural Engineering. 
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plastic and elastic buckling stresses for the assumed half 
wave length L or ratio 8 = L/b. The plastic buckling 
stresses follow from the author’s theory of plastic sta- 
bility. If the coordinate axes coincide with the prin- 
cipal axes of the original state of stress (before buckling), 
the pertinent differential equation is® * 


Ow Ow O*w 
EI| A 2(B 2F) - D— 
| Ox! oT Ox*dy? +r | ° 


oy' 
O*w O*w 
P; Sx? + P, — 0 (36) 


where P, and P, are the compressive forces per unit 
width in the X- and ¥-direction, respectively. 
Furthermore, the strain energy in plastic buckling is* 


l r ow \? Ow O*w 
V, = -El A ) 2B 
2 / J | (= + Ox? dy? . 
ow? o*w \? 
D ( ") + 47 ( = ) | dx dy (37) 
Oy’ Oxoy i 


Values A, B, D, and F are given by Eqs. (21) and (21a) 
of reference 2. They refer to the actual equivalent 
stresses og from Eq. (28) (where in this case 7,, = 
0), computed from the critical loads P,, (rather than 
from Py or P;). Values // in Eqs. (36) and (37) are 
equal to 


EI = (1-—v?)D, or (1 — v?)D, (38) 
respectively, as far as Py or P; are concerned, where D 
and D, are given by Eqs. (20) and (21). 

For Cases A and B (Table 1, Part I), values n follow 
directly from references 3 or 4. For Case B, with 
Poisson’s ratio v = 0.3, 

(A/B?) + 2(B + 2F) + Ds 
ne = 0.91 ~ : (39) 
(1/62) + 2+ @ 
In Case A.1, however, no state of plane strain in Y- 
direction (Table 1) will occur in Case 1, so that for P; 
in Eq. (19) the plastic reduction coefficient /,/F for a 
column obtains, where /, is the tangent modulus. 

For the plastic buckling stresses in Cases C, D, E, and 
F, Eq. (36) yields transcendental buckling conditions‘ 
analogous to those in the elastic range. For the elastic 
range, however, Bleich® gave accurate explicit expres- 


sions for the buckling coefficients ky—namely, 
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Plate subjected to stresses oz, oy, and rzy in its plane, 


Fic. 14. 


ko = (1/8) + p + GB? (40) 


where p and gq are given in Table 1. Hence, in a simi- 

lar way as in reference 6, it may be concluded that, 

analogous to Eq. (39), 7 for Cases C and D may be ex- 

pressed as 

(A/B*) + p(B + 2F) + gDB? 
(1/8) + p + gB? 


nc,p = 0.91 (41) 
The results of this expression are in good agreement with 
the exact values. 

Although it is evident that in Eq. (41) the first and 
last terms of the numerator, which refer to bending 
by \/, and M,, respectively (for X- and Y-directions 
see Table 1), should contain the pertinent coefficients 
A and D, respectively, it is not directly evident that the 
multiplicand of » in the middle term should be (B + 
2F). The plastic buckling stress may be calculated 
approximately by the energy method by equating the 
internal work V; from Eq. (37) to the external work 
V.; hence, from 


(42) 


where, for Cases A-F 


Ve = (1 2)P.S S (ow Ox)* dx dy (43) 


In Cases C and D, where all edges are supported, 


SS (Q’w/dx?) (0°w/dy") dx dy = 
SS (O’w/dxdy)*dx dy 


so that it follows from Eq. (37) that the plastic buckling 
stress and, hence, 7 will indeed contain the term (B + 
2F). In Cases E and F, however, the above terms are 
not equal, and, since torsion plays a major role here 
the plastic buckling stress may contain a larger con- 
tribution by F, related to torsion and a smaller one by 
B, which refers to the lateral contraction. Indeed, for 
Case E, it follows from Eq. (42), letting 
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w= Cy sin (rx/L) 
that 


ne = 0.91 [(A/8?) + 1.21F] /[(1/8?) + 0.425] (44) 


For Case F, assuming 
w = w[l — cos (ry/2b)]| sin (rx/L) 
Eq. (42) yields 


P, = [(A/B%) — 0.3B + 2.2F + 


0.138D B?|r*kI/b? (45) 
which gives for the elastic range, where A = D = 1.1, 


B = 0.33, and F = 0.385, 


kyo = (1/87) + 0.68 + 0.1388" 


instead of the more accurate value 
ko = (1/ B?) + 0.57 + 0.125," 


from Table 1. This is due to the fact that, in general, 
the energy method yields values that are too high. 
Hence, Eq. (45) may be improved somewhat by ad- 
justing it to the more accurate values by multiplying 
the coefficients of B and F by 0.57/0.68 and that of D 
by 0.125/0.138, so that for Case F 


(A/8*) — 0.25B + 1.84F + 0.125D,° 
(1/8?) + 0.57 + 0.1258? 


ne = 0.91 (46) 


In Cases I and J, where P, = P, = P, for Case 2, 
coefficient a and, hence, P2 in Eq. (30) are independent 
of the form in which the plate buckles, so that it buckles 
in that form that makes Py and P; minimum—that is, 
in the same form as a homogeneous plate. The plastic 
buckling stresses for these cases were calculated in refer- 
ence 7, from which the reduction coefficients » = 0.91A, 
as given in Table 1, follow immediately. 

In Fig. 7, values n for Cases A, I, and J are given for 
avional. The dotted curve represents values y for 
Cases Ao, I, and J if in Eqs. (21) and (21a) of reference 
2 Poisson’s ratio v in the elastic range is assumed to be 
equal to 0.5, as is done in Stowell’s formulas. (This ap- 
proximation was first introduced into the present au- 
thor’s theory® by Ilyushin.) It is seen that with this 
approximation the » values for Case A, agree very well 
but that for Cases I and J they are underestimated by 


about 20 per cent. 


In order to derive the plastic reduction coefficients for Cases G, H, and K, the more general case, where the plate 


is subjected to a plane state of stress o,, oy, Tz, (Fig. 14), will be considered first. 


For this case, where the coordi- 


nate axes do not coincide with the principal axes of the original state of stress, the relations between excess stresses 
and excess strains with incipient buckling are given by Eq. (8) of reference 7 and follow also from Egs. (12) of 


reference 2. 


Here, do, = dry. = dtzz = 0, while only de,, de,, and dy,, need to be known, so that 
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de, = (0€,/00,)do, + (O€,/O0,)do, + (O€,/OTz,)d Tz, 
de, = (0¢,/O00,)doa, + (O€,/Oc,)da, + (O¢€,/Or,,)dr,, 
dy zy = (OYzy, Oa,)do, = (Oxy Oa,)do, + (O¥ry a 


The partial derivatives follow from Eqs. (18) of reference 7 or also from Eqs. (13), (14), (17), and (18) of refer- 


ence 2. According to the latter 


G, = (1/3)E,; tan ¢, = (2/3)tan ¢ 


E E EE, EE, (48) 


/ = = > _-_ = 2 oa" tan = s Ee 
ae ae ee et eee ee 


“ft 


) 


in which £, and :, are the secant and tangent moduli. Furthermore, for the considered plane state of stress 


a, = (20, — a,)/3; & = (20, — a;z)/3 
6% = —§¢, + d,) a: a = = (2 3)oR" 
(49) 
where | 
Op = (o,* + oy* — orey + STzy’) 


gx is the equivalent stress according to the shearing energy criterion of plastic deformation as also given in Eq. 
(28). With these data it follows, for example, from Eqs. (13), (14), (17), and (19) of reference 2 and more directly 


from Eqs. (10) and (18) of reference 7 that 


Oc,/Oo, = (1/E,) + [(20, — o,)?/(40%°E,) | 


O¢,/O0, = O¢,/O0, = —[E — (1 — 2v)F,)]/(2EE,) + [(20, — o,) (20, — o,)|/(4082E;) 
Oc,/Oo, = (1/E,) + [(2¢, — o,)*/(40,°E,) ] 
0¢,/OT2 = OYz,/O0r = 37,,(20, — o,)/(20°E,) (30 
O¢,/OTr = OYz,/O0, = 37,7,(20, — o,)/(20n7E,) 
Oy7,/Orr = [BE — (1 — 2W)E,)/(EF,) + [97,,2/(on°E,) | 
where 
E, = E,E,/(E; — £:) 
From Eqs. (47) it follows inversely, denoting differentials by a prime, such as da, by a,’, 
a,’ = Cner’ + Crey’ + Cisyry’ 
oy = Cree’ + Cre,’ + Casyr,’ (5 
Try. = Cer’ + Cre,’ + fia 
With 
€,) = —20°w/Ox"; €,' = —20°w/Oy"; y,,' = —220°w/dxdy (52 
in which zis the distance of the pertinent point from the middle plane of the plate, Eq. (51) yields 
a, = —(Cy20°w/Ox?) — (Cy20°w/Oy?) — (2C\320°w/Oxdy) | 
oy) = —(Cy220°w/Ox*) — (Cxsd*w/dy?) — (2C2320°w/Oxdy) (53 
Try) = —(C320°w/Ox?) — (Cx320°w/Oy?) — (2C3320?w ana 
Combining the above with (Fig. 2) 
*H/2 *H/2 *H/2 
M,’ J o,2d&: M,' = / o,2dz; M,,’ = — / Try 2 a2 
J -H/2 J -H/2 J -H/2 
yields 
M,’ = —I1[(Cy0*w/dx?) + (Cy20°w/dy?) + (2C\30°w a 
M,’ = —I[(C20°w/Ox?) + (Cxd*w/Oy?) 4- (2C230°w/Ixdy) | (54) 


My,’ = +1[(Ci:0°w/Ox?) + (Cx30°w/Oy?) + (2C3302w/Oxdy) } 


(47) 

















tk 


(47) 


refer- 


(48) 


( 19) 


Eq. 
‘ctly 


D0) 


— 


























TABILITY OF SANDWICH PLATES 798 


Since here the unit compressive forces P, and P, are considered as positive, the condition of equilibrium is accord 


ing to reference 9, page 305 


07M,’ 0? M,,’ 07M,’ O7w O*w O7w 
se: 2 + = a es ad agi + rs 
ox* Oxdy Ov" Ox? Oxoy Ov’ 


Insertion of Eqs. (54) in the above gives the pertinent differential equation 


O'w O'w O'w O'w w 
1| Cs Ox! + 43 + 2(Cyw + 2C33) —— -- 40s; - + Cx + 
x 


" Ox8Oy Ox*Ov? Oxdy* oy' J 
P, ‘7 — 2P,, = + P, 4 = (0 (55) 
ox* Oxoy Ov? 


The same Eqs. (51), (54), and (55) and the pertinent values for the coefficients Cy, Ci, etc., may be found by cal- 
culating values A, B, D, and F from Eqs. (21) and (21a) of reference 2, assuming the X’- and Y’-axes in the direc 
tions of the principal stresses. Since stresses and strains transform as the components of a tensor, subsequently, 
the pertinent transformation equations" are used to transform the stress-strain relations to those for X- and VY 
axes parallel to the edges of the plate. 

The strain energy I’, in the buckled plate is (reference 9, page 306) 


V,= -(1/2)S S [(M,'0°w/dx?) + (M,'0°w/dy?) — (2M,,'0°w/Oxdy)] dx dy (56) 


Insertion of Eqs. (54) yields 


/f lc ‘ 4] Pe a (y+ Sere 
> " Vox? * Ox? Oxoy Oxdy ee x? > Oy? 


_ Ow Ow , (ow? a 
4 Coz + Cx ( ) dx dy (57) 
Oxdy Oy" oy" . 


Using oblique coordinates Xo, vy according to Fig. 15, as was done also in reference 11, so that 
x = Xo + Yo SiN g; ¥y = Ww cos ¢Y (58) 


Eq. (57) transforms into 


V, aie? + r Sf} i n = LCha + 2( Cie + 2C33)a? — 4C. 93003 + C. 220 ‘] (cos ¢) (0°w Oxo")? aa 


“41 -_ (Cre + 2C33) a 4. 3Cx30? sat Cra? | (O°w Oxo”) (O72 v ‘Ox pOVo) + 4(C33 = 2Co3a + Cora *) x 
(1/cos ¢) (0°w/Ox0Vo)? + 2(Crw — 2Cx3a + Cra’) (1/cos g) (0°w/OxX_?) (O?w/Ay¥_?) + 
1(Co3 — Coa) (1/cos*g) (0°w/Ox OV) (02w/AV?) + Cxw(1/cos? g) (O°w/ODyo2)?} dxo dyo (59) 
in which* 
a = tan¢ (60) 
Considering Case K of Table 1, where only o, and r,, differ from zero, the external work is (reference 9, page 314) 


Ve. = (1/2) S S |P(Ow/dx)? — 2P,,(0w/dx) (Ow/dy)] dx dy (61) 


which for oblique coordinates transforms into 


~) ; Ow Ow aw \ ? 
-sf f 4p. ) (=) cos ¢ — 2P,, kK > _ (5 ) sin lj dx dyo (62) 
.w 0 Xo Ve ) Xo 


Letting, in Eqs. (59) and (62), 


W = Wy COS (4Xo/L) cos (aryo/bo) 
the condition V; = V, yields, with P, = yP,,, 


‘ae = (r7I/b?) (¥ + 2a) {[Cu - 4C\3a + 2(Ce + 2C23) a? — 4Cx 30° + Cxa4] (1 B?) + 
2(Cyw + 2Cs3 — 6Cox3a + 3C2a?) + CxB?} (63) 


where 8 = L/b. From Eqs. (51) it is clear that in the elastic domain 


* This @ has, of course, nothing to do with « for Case 2 in Eq. (30) and Table 1 








794 JOURNAL OF THE AERONAUTICAL SCIENCES DECEMBER, 1951 
Cu = Cy = E/(1 — v?); Cy = vE/(1 — v?); | 
" , ’ . (6 
C33 = E/2(1 + »); Cy = Cx = 0 ( 4 


Hence, it follows from Eq. (63) that for a fixed value of 8, as is considered in case of sandwich plates, the reduc- 


tion coefficient for plastic buckling for Cases 0 and | is, with vy = 0.5, 


O.91[(y + 2a)—'} [Cu 


le \E(y + Pa) (1 


1Ch;a -+- 2( Cre + 2C%33) a? 


+ a”)? 


— 4Cx3a* + Cra*|B ules 
2(Ce + 2C33 = BCo3a + 3Cxa?) + Cx» 8?! Lest 


( 65) 


B+ + 2(1 + 3a") + B13 win 


where the subscript (min.) indicates that a = tan ¢ has to be chosen such that both numerator and denominator 


obtain minimum values. 


Considering the case that a, and 7,, contribute equal 
portions to the equivalent stress o, in Eqs. (49) and, 
hence, that ¢. = t4V 3 0ry = P./P. = Ge) Ty 


V 3, it follows from Eqs. (50), since o, = 0, that 


Oc,/Oo, = (1/E,) + [1/(2E,)] 
Oc,/Oo, = Oe,/O0, = —j[K — (1 2v)F,| 

(2EF,)| [1 /(4F,)] 
Oc,/Oo, = (1/F,) + [1/(8E,)] 
O¢,/Orr, = OYz,/007 V 3/(2E,) (66) 
O¢,/Orry = OFz,/O00, —V3/(4E,) 
Oy1/Ory = {[BE — (1 — 2)E,|/(EE,)} 4 

[3/(2E,) | 


These values were evaluated for the case that the equiv 
from Eqs. (49) reached 
0.1675, 


alent stress og = (6,7 + 37,;,7) 
a point of the stress-strain diagram where e 
tan g = 0.294E or FE, = 0.857E and EF, = 0.227E. 


From this it follows from Eqs. (47) and (51) that 


Cy = O.700E; Crp 0.3694 
Cs = —O.1S1E; Cu 0.950E (67) 
Co; = 0.037E; Css 0.1964 


With these values Eq. (65) yields, for an assumed value 
8 = L. 


Ne = 0.744 (GS) 
From Eq. (65) the reduction factors for Case G, pure 


shear, follow by setting y | ak ke 0. Further 


more, from Eqs. (50), since o, Oy 0 and o,” 
Sieg’: 

Oe,/Oa, O€,/ Oo, I/E 

O€,/ Oa, Oe, Oo, (1 2v)/ (2E)] 


| | one (69) 
E| + (3/E, 


Hence, from Eqs. 


OV x5/ OF ay (1 2v) 


while the other derivatives vanish. 


(47) and (51), 


Cu = Cos = A E> Cie = BE; Co = FE (70 


where A, B, and F are the coefficients given for Case G 
in Part I of this paper (JOURNAL OF THE AERONAUTICAL 
SCIENCES, May, 1951) which were derived directly for 
this case in references 3 and 7. The other coefficients 
C are zero. 
ing Eqs. (51) with the first two and the last one of Eqs. 


The relations (70) follow also by compar- 
(21) of reference 2. Hence, for Case G, Eq. (65) yields 
the expression for y,, which was given in Eq. (26) of Part 
I. In Fig. 13, ng is given for avional. Stowell found in 
reference 21 that for homogeneous plates the 7 values for 
Cases G and H are practically equal. Since his for- 
mulas are identical with ours, except that he sets y = 
0.5, it may be concluded that this identity will remain 
true if vy has its actual value of about 0.3. Since, more 
over, n; varies little with 8 (Fig. 13), it may be as- 
sumed that, also for a definite value 8, the reduction 
coefficient ny is equal to ng. 

Since n; 1S never appreciably larger than its value 
for B 
curately simplified by assuming 6 = 0. 


0 (Fig. 13), Eq. (26) may be sufficiently ac 
This yields 
A(l +- at) + 


2(B + 2F)a? 


NG NH (0.296) a) 


where 


(26a) 

















YV Ny, 


[rausformation from coordinate axes Y, 
nate axes Yo, | 


}° to coordi 


Fic. 15 
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For the same equivalent stress og as used above in 
calculating 4, and, consequently, for 7,, on/V 3, 
values A, B, and F as given for Cases G and H in Part I 
are (compare also references 3 and 7) 

A= 09600: B = 0.316; F = 0:078 
so that from Eq. (26), with 6 = | 


as = 0.701 (71) 
The simplified Eq. (26a) yields ng = 0.698. 

ne (for Case B) is also 
For this 


In order to check Eq. (27), 
calculated for the same equivalent stress. 
case of pure compression, value y in Eq. (65) is infinite, 
and Cy, = AE, Cy» = DE, Cy = BE, and C3; = FE, 
Hence, Eq. (65) 
For the same 


while the other coefficients vanish. 
reduces to nz from Table | and Eq. (39). 
equivalent stress o,, hence for ¢, = op, here [reference 
7, Eq. (27)| A = 0.421, B = 0.426, D = 0.938, and F = 
0.322, so that Eq. (39) yields 
Since the state of stress in Case K is between those in 
Cases B and G, it should be possible to express nx in 
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nx and ng. A plausible assumption is to allow for nz 
and ng an importance proportional to the contribution 
of o, and 7,, to the shear energy or to og”, so that from 


the last of Eqs. (49) 


Ce" 4 Se 
ig renee a 
8 a2? + Bry? Or° + Stry° 
or, with o, = VTews 
nk = (¥*np + 3ng)/(y? + 3) (27) 


where all » values refer to the same equivalent stress 


on and to the same ratio 6 = L/b.* 


From Eqs. (68), (71), and (72) for the given og and 


8 


nx = 0.744; ng = 0.701; ng = 0.796 
With the latter two values Eq. (27) yields 
nx = (1/2) (np + neg) = 0.748 


so that the error is '/» per cent. 


PLASTIC REDUCTION COEFFICIENTS FOR HOMOGENEOUS PLATES 


For homogeneous plates, where 8 adjusts itself such as to minimize the buckling stress, Eqs. (41), (46), (26), and 


(65) lead to the following reduction coefficients: 


ne.n = 0.91[2(gAD)'? + p(B + 2F)]/(2Vq + P) 


where ~ and g are given in Table 1. 


ne = 0.504(V AD + 2.6F — 0.35B) 


nc = nn = 0.161[a-"({A[A(1 + at) + 2 


0.91 [(¥ + 2a) . (} Co[Cu = 4Cx3a + 2( Cie + 2C33) a? = 4C3303 + Cora] }' 7 + Cie zi 


= > 


(73) 

(74) 

(B + 2F)a?)}'? + 34a? + B+ 2F)) min (75) 
2C33 ~~ 6Co3a + 3Cx2a”) |min { _- 

(7460) 


{2E[(y* + 2)" — 4]} 


In Eqs. (75) and (76), a has to be chosen such as to make the pertinent expressions a minimum. 


In order to check Eq. (27) for this case, n, and ng were 
calculated from Eqs. (75) and (76) and ng from Table 
5 of reference 2, 


nn = 0.455(V AD + B + 2F) 


for the same equivalent stress o, as used above and 


with y = V 3in Eq. (76). The results aret 


ne = 0.744; ng = 0.715; np = 0.773 


tT In reference 7, page 33, the exact value of ng for the same 
point of the stress-strain diagram was found from the differential 


equation of the plate. It is 3.88/4.82 = 0.701 


With the latter two values Eq. (27) gives 


nx = (1/2) (ne + ne) = 0.744 


This is in exact accordance with value nx from Eq. 

(76), so that Eq. (27) can be used instead of Eq. (76). 
Formulas (73), (74), (75), and (27) form an extension 

of Table 5 of reference 2. The formulas for Cases C 


* About at the same time that the above derivations were pub 
lished,! the case of plastic buckling of a homogeneous plate under 
combined shear and compression was considered by Stowell," 


who came to another way of approximation 
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and F differ somewhat from those of Table 5 because 
they were adapted to Bleich’s values p and q. 


THE BUCKLING COEFFICIENT (kj), FOR CASES 0) AND | 


For the elastic range, with C values according to 
Eq. (64), Eq. (65) for Case K transforms to 


Fas = [r*EI b7(1 = v) |(Ro)x ii 


[r?EI/b2(1 — v*)| (y + 2a)! X 
[((1 + a?)*B-? + 2 (1 + 3a?) + B?] (77) 
where a = tan ¢ has to be chosen such as to make P,, 


minimum. With y = P,/P, = V3and B = L/b 
1, 0.5, and 0, respectively, this equation yields values 


(Ro) = 2.02; 3.20, 0.51/ 8? (78) 


respectively. On the other hand with y = 0, Eq. (77) 
gives the buckling coefficient (ko),; for Case G, yielding 
for the same £6 ratios 


(Rog = 5.83, 9.85, 1.54/B? (79) 
while, since (ko)g = (8~! + £B)*, for B = 1, 0.5, and 0 
(ko) = 4, 6.25, l B? (SO) 


Inserting Eqs. (79) and (80) in Eq. (25), one finds 
(ko) « = 2.04, 3.25, and 0.51/B?. 

This is in good accordance with Eq. (78). Hence, 
the accuracy of Eq. (25) for (ko), is satisfactory. 


THE BUCKLING COEFFICIENTS @ OF CASE 2 


A rather complete derivation for all cases in Table | 
was given in reference 1. However, the derivations for 
Cases A,'? B, C, D, I, G, H,'* and K' were already 
given in the author’s earlier papers, so that only Cases 
E, F, and J will be considered here. 

The differential equation of the plate for Case 2 is! 7 


(¢ + h)? G. (= =~) _ Pp, Ow _ on , 
t Ox? ov? Ox? " Oy? 
ZF xe : ae 0 (S81) 
Oxoy 


where G, is the pertinent modulus of rigidity of the 


core. With the coordinate axes as in Table 1 for 
Cases E and F, the unit forces P, and P,, are zero. 
With 
w= Y sin (xx/L) (S82 
in which Y is a function of y alone, Eq. (S1) yields 
VY" + (x?/L*) {tP,/[(¢ + h)*G.] — 1}¥ =0 (83) 
where a prime denotes differentiation with respect to y 
The solution is 
Vv Ci sin gy + Cy cos gy S4) 
where 


¢g = +(n/L) }tP,/[(t + h)°*G, 1} (85) 
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At the supported edge (y = 0) for both Cases E and F 
the boundary condition is 


(Y),.0 = 0 (86) 


while at the free edge (y = 5) the shear forces Q, vanish 
Assuming that at the short edges (parallel to Y-axis, 
Table 1) the faces are prevented from a displacement in 


the Y-direction, from Eq. (7) of reference 13 
Q, = [(t + h)? G,/t] (Ow/dy) 
Using Eq. (S82), this yields 
(Y’).., = © (87) 
From Eqs. (84), (86), and (87) 
C. = 0; cos ¢b 0 
Hence, for a minimum value of P,, from Eq. (85) 


P, = [1 + (L?/4b?)] (¢ + h)°G./t = 
[1 + (B74) |(¢ + h)°G,t 


so that from Eq. (30) for Cases E and F 
Qe = ar = 1 + (87/4) (SS 


For Cases I and J, where P,, 0 and P, = P 


P, Eq. (S81) transforms to 


i [(t + h)°G,/t] — Py [(0?w/dx?) + 
(0?w/dy?)| = 0 (89) 


whence 
P, = P, = ((t¢ + h)?/t\G. 


independent of the form of the deflection surface, and, 
from Eq. (30), 


= ] (90 


Since here, consequently, P2 in Eq. (11) is constant, for 
these cases a sandwich plate buckles in such a form 
that Py) and P; are minimum and, hence, in the same 


form as a homogeneous plate. 
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Aerodynamic Coefficients of an Oscillating 
Airfoil in Two-Dimensional Subsonic Flow’ 


R. TIMMAN,? A. I. VAN DE VOOREN,! AND J. H. GREIDANUS** 
National Aeronautical Research Institute, Holland 


SUMMARY 


A method yielding a direct solution in terms of known func 
tions is presented for the aerodynamic forces acting on an oscil 
lating wing in a two-dimensional subsonic flow. The solution is 
based on the linearized equations of motion. Lift and pitching 
moment coefficients are calculated for a translation and for a ro 
tation of the airfoil. Results are presented for five different 
values of the Mach Number (0.35, 0.5, 0.6, 0.7, and 0.8) and for 


values of w ranging from 0.1 to 3. 


List OF SYMBOLS 


rectangular coordinates in the physical plane 


Yrv= 

l = semichord 

1 = air speed 

( = speed of sound 

g = Mach Number, v/c 

t = time 

» = frequency 

w = reduced frequency, vl/v 

¢ = acceleration potential 

& = velocity potential 

po = density of the undisturbed airflow 

A = amplitude of translation (positive if downward ) 

B = amplitude of rotation about mid-chord (positive if 
trailing edge downward ) 

K = total force on airfoil (positive if downward ) 

M = total moment about mid-chord point (positive if tail- 


heavy ) 
Subscripts denote partial differentiation with respect to the 
subscript variable. 
A bar above a letter denotes an amplitude. 


INTRODUCTION 


— VALUES OF THE aerodynamic coeffi- 
cients of an oscillating wing in two-dimensional 
subsonic compressible flow have been calculated, for the 
first time, by Possio.'. These calculations were based 
on his own integral equation for the chordwise pres- 
sure distribution. They have been repeated and ex- 
tended by Frazer.** No analytical solution being 
known at that time and the kernel of the equation 
being highly complicated, the rather crude and pro- 
visional direct numerical methods used are not apt to 
accuracy. More careful and elaborate 
numerical results from 


yield high 


methods to obtain Possio’s 


equation, still without solving it analytically, have been 


* This investigation was sponsored by the Netherlands Air- 
craft Development Board 

+t Research Mathematician, 
matics at Mathematical Centre. 


Chief Section Applied Mathe- 
t Research Engineer 
** Chief Section Flutter and Theoretical Aerodynamics. 
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developed and employed in Germany by Schade‘ and 
by Dietze.2 The method of Schade re- 
duction of the integral equation to a finite set of alge- 
braic equations, and the method of Dietze is an itera- 
The maintenance of satisfactory ac- 


involves a 


tion method. 
curacy should, however, still be considered difficult in 
view of the fact that no rigorous limits to errors can be 
indicated. 

During the war, Timman® developed a rigorous ana- 
lytical method to solve the equivalent boundary value 
problem to the wave equation based on expansions in 
terms of Mathieu and modified Mathieu functions. 
By this method, complete analytical solutions applying 
to a plain wing and to a wing with an aerodynamically 
balanced control surface have been elaborated in col- 
laboration with van de Vooren.’ This work is now 
being used to obtain new numerical results. The ex- 
tensive computations required are being performed in 
the Computation Department of the Mathematical 
Centre at Amsterdam under the direction of van 
Wijngaarden and with the assistance of Scheen and 
Berghuis. 

The results of the first part, applying to the plain 
wing, have just been completed and are given in this 
Other results will follow as soon as they are 
ready. All tabulated values should be accurate up to 
one unit in the last decimal supplied. 

Graphs and tables will be preceded by a short review 
of the mathematical principles involved. Other ap- 
proaches to the subject involving Mathieu functions 
have been given by Reissner and Sherman,'* Haskind,"” 
Billington,’ and Reissner.'* 


paper. 


AERODYNAMIC THEORY’ 


The Mathematical Problem 
It is well known’ that the linearized equation for 


the acceleration potential 
¢(x, ¥, f) = g(x, ye” 


of an oscillating airfoil in a two-dimensional subsonic 
flow can be reduced to the time-free wave equation 


(Pockel’s equation) 
Uxx + uUyy + y/o) GB") lu 0 é 
involving independent variables 


X = x/V 1 — 8B: 
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(Prandtl-Glauert transformation), and the auxiliary 

function 

u(X, Y) = gXV 1 — B’, Y) exp (—wBX/cV 1— B?) (3) 
The problem consists of finding a solution to Eq. (1) 

which: (a) satisfies the boundary condition 


Lo : , x—€& 
> go(—, ¥) exp iv(t — ) dé = w(x, t) 
0 OY J —. 


(4) 


on the segment —/ < x + / of the axis of x, the verti- 
cal velocity w(x, f) following from the assumed motion 


= vl 
v(x, t) = H(n)e’ 


of the points of the centerline of the airfoil by substan- 
tial differentiation 


w(x, t) = v(x, t) + vy,(x, t) (5) 


(6) is regular everywhere in the closed region of the 


x, y plane, cut open along the segment —/ <x < + / 


of the axis of x, except in the point x = —/, y = O rep- 
resenting the leading edge of the airfoil, where it should 
become infinite like R-* (R = distance to leading 


edge) and except at points where the (piecewise con 
tinuous) function w(.v, /) is discontinuous. * 

This solution consists of a regulart part uniquely 
determined by the condition 


(Prex Vy UW, t Os 4 y = 0) (6) 
following from Eq. (4) by substantial differentiation, 
and a singular part uniquely determined, safe for a 
multiplicative constant, (1) by the condition 


( l<xSQ4l 


p ly QO 
(¢ sing./ 4 { y sO 


(2) by the requirement stated under (b) about the 
nature of the singularity in the point x b ¥ 0. 

The multiplicative constant is determined by the 
requirement that the complete solution should satisfy 
Eq. (4). 

The regular solution, being finite at the leading edge 
(which implies that pressure and velocity are finite 
there), corresponds with shock-free entrance of the 
flow. Both parts of the solution imply vortex sheets 
of the axis of x, the evolution 


} 


along the part x > 4 
of which is implicitly determined by the imposed con 


tinuity of the pressure field. 


The Regular Solution 


Introducing elliptic c ordinates &, n, 


* Such points, representing breaks in the airfoil’s centerline, 
do not occur in the case of a plain rigid airfoil 

+ A solution is defined to be regular when it is continuous at 
the boundary, continuous together with its second derivatives at 
all inner points of the region, and regular at infinity; compare 


Kellogg.’ 


x : l 
=X = cosh & cos 7 
V1—- 6 V1— 8 
‘ y ' , 
y= Y= sinh € sin 7 
V1l- s# 
Os &€ O n < 29 (S 


adapted to the boundary value problem under con 
sideration, the wave equation changes into 


Ure + U,, + 2k*(cosh 2 — cos 2n)u = 0 (9 


k = (1/2)6Q, Q = w/1 — B 10) 


and the boundary condition (6) for the regular part of 


the function u becomes 


u(O, 4) = — @, (O, n) — tw®(0, n) sin n | x 
Vi1— p 
ex — 1B°2 cos n (1] 
with 
w(t, n) = WE, noe 
Putting 
u(é, n) = Uy1(E)u2(n 


and separating variables, the modified Mathieu and 
Mathieu equations 


(Ui)ze — (a — 2k? cosh 2&)u, = O 12 


(U2) + (a — 2k? cos 2n)u. = O 13 


are obtained. 

Since in the variable » the acceleration potential 
(pressure) is odd and periodic with period 27, the rele 
vant solutions of Eq. (13) are the functions, usually" 
denoted by se,(n). The subscript » determines the 
order of the corresponding characteristic value a = ),,. 

The relevant solutions of Eq. (12) are the functions 
Ne,‘ (&), since these solutions represent waves moving 
away from the origin with a speed c. 

Hence, the regular solution of Eq. (1) can be written 


as 


2. anNe,(&)sen(n 14 


where a factor v* has been added to make the quantities 
a, dimensionless. 

By developing the function “,(0, 7) from Eq. (11 
in a uniformly convergent series of Mathieu functions 
sé,(n), the coefficients a, can be determined by com 
parison with the function u,(0, 7) following from Eg 


(14). The regular solution is thus completely known 


The Singular Soiution 


The singular solution can conveniently be derived 


from the function 


(g, n = sin n/(cosh € + cos n 
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representing the acceleration potential of the steady 
compressible flow about a fixed flat airfoil at some non 
vanishing incidence and thus showing a singularity of 
the required type at the leading edge (& = 0, 7 = 7m). 
Multiplying by e’”, the function does not remain an 
acceleration potential, because it does not satisfy the 
wave equation [Eq. (1)] replacing Laplace’s equation 
in the case of unsteady flow. The function is modi 
fied into a solution of the wave equation by a series of 
correction terms, preserving the condition of vanishing 
normal derivative at the contour — = 0. 
Using the series expansion, valid for — > 0, 

se = 2 >> (-1) "e "sin nn 
cosh £ + cos 7 care 


the modified function is formally written as 
j sin 7 


lcosh — + cos 7 


oS Hu ro(2) (t eo: { . 
- »® c=—a7 a, Ne‘), (é)se,(n) — e ™ sin nn i (15) 
n=1 
which for — > 0 is a solution of the wave equation. The 
coefficients a, are determined by the condition 


u:(O, n) = O. 
For the purpose of determining a@,, the expansion 


sin mn = >, B,” sen(n) 
l 


m 


is introduced, and the order of summation is inter- 


Eq. (15) then formally changes into 


u(é, 9) = v - —2 >> (-1)"X 
lcosh € + cos 7 n=1 


changed. 


i sin 7 


a,Ne®,(é) — >> B,'” e~™ S€n(n) ¢ 

m=1 
and the coefficients a, are determined by imposing the 
condition u;(0, 7) = 0 on each of the terms in this se- 
ries. The resulting expression becomes 


u(é,) = v +2 >> (-1)"X 


(cosh & + cos 7 Pigg 


j sill 


Ve,,°?)(&) sa r 
se,’ (0) > ne 2 


—~ 


{ 
+ sé,(n) (16) 
Ne,, (2)(() m=1 . 


The series of correction terms can be shown’ to be 
uniformly convergent, together with its first deriv- 
OU s 


atives in the region & n < 2n. 


The Resulting Potential 


Writing 
u(t, 7) = Urer (, 9) + Gotlsing (E, 2) 


R.H.S. Eq. (14) + do R.H.S. Eq. (16) 
the arbitrary constant dy) is to be determined by going 
back to the velocities at the boundary and imposing 
Eq. (4) in one arbitrary point —/ S x tl,y = 0. 


In this calculation, the singular function (16) again 
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causes a slight difficulty in that, for this function, inter- 
change of differentiation and integration in Eq. (4) is 
not admissible. This difficulty can be eliminated by 
aid of a partial integration (see reference 7). 

The function u(£, 7) now being known, the accelera- 
tion potential and the pressure distribution can be at 
once written down, and aerodynamic forces and mo- 
ments can be obtained by the appropriate integrations. 


NUMERICAL CALCULATIONS 


The complex dimensionless derivatives k,, k,, ,, Mp, 


defined by 


K = mpolv*e’" (Ak, + Bh) 


M = mpyl?*v7e'"(Am, + Bm,) 


were calculated for five values of B—viz., 8 = 0.35, 
0.5, 0.6, 0.7, and 0.8—while, for each 8, ten different 
values of w, ranging from w = 0.1 tow = 3, were taken. 
The values of w were chosen in such a way that, by Eq. 
(10), they led to a limited number of simple k-values, 
since k is the parameter determining the Mathieu 
functions. 

The functions se,(y) have been obtained by calcula- 
tion of the coefficients B,,””’ from the Fourier series 
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Fic. 3. Aerodynamic coefficients as functions of the reduced 
frequency w with Mach Number £8 as parameter. Points calcu- 
lated for w < 1 have not been marked separately. 


frequency w with Mach Number £ as parameter. Points cal 
culated for w < 1 have not been marked separately. 


sen(n) = >. B,'” sin nn (17) 
-n=1 

For small values of g(= k?*), the coefficients B,°””’ and 
the characteristic values }, were calculated from their 
power series!! to g. For g 2 0.3025 these quantities 
were calculated from their expansions in factorial 
series, |! which form an application of the general theory 

of factorial series as developed by No6rlund.!” 
The modified Mathieu functions Ne, ?(£) are written 


as a linear combination 
Ne, (&) = Se,(~) + aGe,(&) (18) 


of the two independent solutions Se,(é) and Ge,(&) of 
Eq. (9). In analogy to Eq. (17), the expansion for 


Se,,(&) is 


Se,(é) = }> B,”? sinh mé 


m=1 


while the following relation” holds for Ge,,(&): 


Ge,(é) = —é&Se,(~) + > gm" cosh mé 
m=0 


(n) 


The coefficients g,,\"’ can, similar to the B,,’"’, be 
calculated either from their recurrence relations or from 
their factorial series. Finally, the factor a, by which 
Se,(~) and Ge,(é) are linearly combined to yield the 
function Ne, )(£), is determined by the condition 
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that the required combination behaves asymptotically 
as H,?(ke*), where //,,” is the Hankel function of the 
second kind. 

When calculating with aid of Eq. (18) the vertical 
velocity at the airfoil, there appear integrals of the type 


—iQ2_ cosh Er * 
f e 122 co 'S6(£) dé 
0 


with /() = cosh mé, sinh mé or & sinh mé. 
that the integrals with f(£) = cosh mé are proportional to 
Hankel functions of the second kind, while the integrals 
with f(é) = sinh mé have an elementary character. By 


It is known! 


partial integration, the integrals with f(&) = & sinh mg 


Approximations 


To obtain additional information about the derivatives in the region of small w-values (i.e., w 


following expansions were established :'4 


© 
— p (2442 +1+2In- — 672+ F(8)|t — (22 — r2)i) 
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can again be reduced to Hankel functions of the second 
kind. 

The coefficients a, in the expansion (14) of the regu- 
lar part of the function u can also be calculated, and it 
follows without difficulties that they are expressible 
in Mathieu functions of the argument 6 = cos~! 8. 

The resulting values for the aerodynamic derivatives 
ka, Roy May My are given in Table 1, while Figs. 1-4 
illustrate the general character of these quantities as 
dependent of 8 and w. In the figures, only points with 
w 2 1 have been marked. The valuesin Table | are 
guaranteed to be accurate up to a few units in the last 


digit. 


< 0.1), the 


k= V1 
Q l d 
m=—-Vi1- a*(24} + In= — = B[2 + F(B)]e — 5 22 — O*)e 
ee € Be rn, : 
ky = —2 + rk + OY — —— § + 28‘ + - # In— + (1 — 28") X 
V1 — B 2 4 2 2 


Q 
(, + In ;) 


{ 


3 { 
= at)F(A) +2 E + In 5 


afr 


2 Boo, 5 
— 5 B’F(B) ‘io 
- (19) 


Q l Q 
(2 dy + 2y + 21n = as B13 + F(B)I¢ = 70") - + 2y + 21n oe B7[1 + rit) | 


l I ifr? l 9 
my, = l ——-@ + aXk— +--+ 
V1 — B 4 4° Ss i 


(, 
3? 
9 


j Q l .) 
(2) y + In = ae B2 (2+ F(B)] 


where y denotes Euler’s constant (y = 0.5772157), 


while F( 8) is given by 
—— vi =e, tui te) 
1 
B 


F(g) = -In- + = 1 
: i B- 


The error in the approximations (19) is of order 
2% In Q. It is estimated that the neglected terms affect 
only the third decimal! digit, provided 2 is smaller than 


0.05. 


Pe ee > Q 
ere = sar y¥+in; + 


Q l ‘ 
- 3) F (B) — E ta a(8) | + 


J 


mm 2 P&S F(a) | 
~ey tremens at (8) () || 


The diagrams show that the derivatives assume an 
oscillatory character if both w and £ are large (k being 
the decisive parameter). This makes interpolation 
rather uncertain when w exceeds the value 1. A further 
investigation of the derivatives in this region will form 


the subject of a separate publication. 


Comparison with Other Results 
Until now, numerical values for the derivatives have 


been given by Schade’ [w = 0, (0.2), 1, 6 = 0.5, (0.1), 
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TABLE 1 


Aerodynamic Coefficients for 


| 


w ee’ Og Ma,’ 
B = 0.35 0.12536 —0.038045 — 0.20720 0.028172 
0.25071 —0.046232 — 0.35499 0.058948 
0.37607 —0.009959 — 0.48879 0.084936 


0.50143 +0 .06580 — 0.62253 0.10905 


0.62679 0.17760 — 0.76220 0.13361 
0.75214 0.32346 — 0.91167 0.16032 





1.00286 0.71222 — 1.253857 0.22520 
1.25357 1.22655 — 1.67687 0.31275 
1.88036 2.94154 — 3.31464 0.68339 
3.25929 5.69807 — 10.6691 2.26557 
B=0.5 0.15 —(0.054862 — 0.24624 0.042830 
0.3 —0.058224 — 0.41738 0.087937 
0.45 —0.003371 — 0.58385 0.13209 
0.6 +0 .09903 — 0.76466 0.18199 





0.75 0.24182 — 0.97049 0.24237 
0.9 0.41792 — 1.21067 0.31685 
1.125 0.72441 — 1.65357 0.46004 
1.650 1.39131 — 3.12396 0.92502 
3 1.50438 — 7.08452 1.48907 
B= 0.6 0.10667 —().05034 — 0.19603 0.034663 
0.21333 —(0.08031 — 0.32744 0.075165 
0.32 —().07722 — 0.44631 0.11429 


0.15588 
0.20283 


—(.04808 — (0.56767 
0.69827 


0.42667 
0.53333  +0.00180 —- 


0.64 0.06767 — (0.84310 0.25724 
0.8 0.18783 — 1.09168 0.35411 
0.96 0.31597 — 1.38548 0.46874 
1.17333 0.46521 — 1.84129 0.63545 
1.70667 0.60534 — 3.05547 0.92597 
3.2 +0 .96358 — 6.33758 1.00499 


B =0.7 0.10929 —0.06566 — 0.20758 0.045313 
0.21857 —0.10485 — 0.33866 0.096150 
0.36429 —0.10752 — (0.49998 0.16528 
0.54643 —(0.06375 — 0.72084 0.26618 


0.33432 
0.42586 
0.50418 
0.57313 
0.645438 
0.64424 


—0.02552 — 0.87102 
— 1.09191 
— 1.32783 
0.10495 — 1.68198 
0.25388 — 2.59601 
0.97723 — 4.64536 


—0.09107 _ 


0.65571 
0.80143 
0.94714 
1.16571 
1.74857 
2.91429 





0.063153 


8=0.8 0.1125 0.21913 
0.225 —0. 14333 — (0.34875 0.12729 
0.3375 —0.16308 — (0.45833 0.18918 
0.495 —(). 16302 — 0.62636 0.26882 
0.585 —(Q.15701 — ().72849 0.30309 
0.72 —0.14703 — ().88543 0.33397 
0.9 —0.1380138 — 1.09148 0.36225 
1.08 —().09809 — 1.28664 0.41640 
1.8 0.46458 — 2.25957 0.45568 
2.4 3.45168 1.23055 


a —().67692 — 


0.8], by Dietze* [w = 0, (0.02), 0.1, (0.1), 0.7, 8 = 0.5, 
0.6, 0.7], and by Turner and Rabinowitz" [w = 0, 
(0.02), 0.1, (0.1), 0.7, 8 = 0.7], the last authors using 
the same method as Dietze. Comparison with these 
values shows that the values of Schade and Dietze agree 
mutually better than they do with the values presented 
here. The differences between Dietze and Schade are 
less than 1'/: per cent in modulus, while the new values 
show differences up to 5 per cent, with, in general, 
greater differences in the arguments. The reason for 
these discrepancies is not known. All steps of the 
present analysis, both in the establishment of the for- 
mulas as well as in the numerical part, were carefully 
and frequently checked, and the authors believe that 
no error has slipped through. 

The greatest differences occur in k, and m, when 
w = | and in k, and m, when w is about 0.5, both for 


an Oscillating Wing Without Control Surface in a Subsonic, Compressible Flow 


ta” k,’ ky,” mM,’ mm," 
0.10315 — 1.68089 +0. 20041 0.83907 —(0. 24038 
0.17540 —1.47408 +0.00944 0.73686 —(). 28289 
0.23925 — 1.38264 —().21186 0.69524 —().31115 
0.30117 — 1.34665 —0.43091 0.68350 —().34275 
0.36340 —1.34317 —0.64462 0.68934 —().38003 
0.42683 — 1.36265 —0.85400 0.70718 —(), 42308 
0.55843 — 1.45668 — 1.26494 0.76829 —(). 52552 
0.69590 —1.61994 —1.66819 0.85423 —0.64961 
1.03275 2.36967 —2.59076 1.13288 — 1.06234 
0.90273 —5.31768 | —2.72070 1.65624 —2.37072 
0.12136 — 1.68380 +0 .24477 0.83345 —0.31277 
0.20082 — 1.47586 —0.00545 0.72258 —().37209 
0.27225 — 1.42102 —(). 26263 0.68800 —0.43239 
0.34195 — 1.44023 —0.50450 0.68552 —0.50786 
0.410385 —1.50998 —().73080 0.69855 —(0.59940 
0.47559 — 1.62300 —(0.93994 0.71826 —(). 70676 
0.55918 — 1.86798 —1.21015 0.74633 —0.89558 
0.58912 —2.69906 —1.48356 0.71801 —1.41513 

—0.48300 —3. 69380 —0.17490 0.26512 —2.19422 
+0 .096452 — 1.87186 +0 .37583 0.92318 —(0.34456 
0.15676 —1.60711 +(). 22084 0.77647 —().41251 
0.20660 — 1.50147 +().03649 0.70842 —().46619 
0.25206 — 1.47193 —0.13766 0.67389 —(0.52852 
0.29418 — 1.48870 —(0.29712 0.65439 —(), 60239 
0.33233 — 1.54088 —(). 48990 0.64064 —(). 68791 
0.37673 — 1.66589 —0.62281 0.61552 —(). 83276 
0.39854 — 1.83983 —(0).75393 0.57364 —0.99169 
0.37720 —2.11168 —(0.82932 0.47683 — 1.20183 
0.08896 —2.63026 —0.57421 0.10857 — 1.52439 
—0.444385 —2.72741 +0) .09652 0.15806 —1.72177 
0.10047 — 1.94357 +0.50111 0.94333 —(), 44332 
0.15523 — 1.63907 +(0).32599 0.75805 —(). 52526 
0.20793 —1.51768 +0 .08661 0.64075 —().63102 
0.24898 —1.54169 —(). 14542 0.53483 —().79581 
0.25586 —1.60498 —(). 24453 0.46137 —(), 90187 
0.23881 —1.71801 —(). 32436 0.33959 —1.03115 
0.193338 — 1.83523 —(0.34523 0.19735 —1.12448 
0.10040 —1.96611 —().30050 —0.00355 — 1.17553 
—(0.02899 —2.07784 —(),. 15406 —(0. 15024 —1.14710 
—(). 72862 —1.98196 +0) .04479 —(). 20001 — 1.52697 
0.10152 —2.00751 +0 .70906 0.93271 —(). 60983 
0.13996 — 1.638890 +(). 49597 0.66763 —(). 70161 
0.15577 — 1.51321 +0 .31454 0.49838 —0.79352 
0.14488 — 1.48420 +0.13713 0.27487 —().90157 
0.12534 —1.50212 +0 .07412 0.14915 —(), 92865 
0.094382 — 1.53969 +0 .02445 —().00200 —().91614 
0.07570 —1.57055 +0 .00753 —().09255 —().87376 
0.06448 —1.57102 +0 .00624 —(0.11260 —(). 89342 
—0.22942 —1.56060 —(0. 28730 —0.35909 —(). 97449 
—0.097835 — 1.43405 +0 .78844 —0.01237 — 1.70254 


8 = 0.8. In the last case (however for 8 = 0.7, since 
Dietze gives no values for 6 = 0.8), the difference be- 
tween Schade and Dietze is also relatively great. 
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A Shock Tube Method of Generating 


Hypersonic Flows 


ABRAHAM HERTZBERG* 
Cornell Aeronautical Laboratory, Inc. 


* SUMMARY 


A method of circumventing certain Mach Number limitations 
in shock tubes is indicated. A diverging approach to the test 
section permits the establishment there of steady flow Mach 
Numbers far exceeding those that could be produced in a constant 
area tube. The advantages of this method of generating hyper- 
sonic flows, with respect to the test section temperature, are dis- 
cussed, and a series of exploratory experiments in a small shock 


tube are described 
DISCUSSION 


— stupiEs! have demonstrated that hyper- 
sonic flows can be generated in shock tubes of 
special design in a manner that eliminates the tempera- 
ture problem* * encountered in hypersonic blow-down 
tunnels. 

The conventional shock tube*~’ consists of a tube 
partitioned by a thin diaphragm into regions of high and 
low pressure. The high-pressure chamber and the Jow- 
pressure chamber are usually tubes of constant area. 
The shock tube is set into operation by destruction of 
the diaphragm. This allows the high-pressure gas to 
expand into the low-pressure chamber creating a shock 
wave that propagates into the low-pressure chamber 
compressing and accelerating the air between the wave 
and the interface (the plane of contact between the gas 
from the high-pressure chamber and the gas of the low- 
pressure chamber). At any point in the low-pressure 
chamber, there is a period of steady flow between the 
times of arrival of the shock wave and of the interface 
in which testing may be carried out. 

Supersonic flows may be created in this region be- 
tween shock and interface by employing strong shock 
Waves (i.e., a large pressure ratio across the shock 
wave), but the test Mach Number is limited by the 
Mach Number that may be created immediately behind 
anormal shock. This Mach Number, which is only a 
function of the specific heat ratio (y) of the gas in 
the low-pressure chamber, approaches a limit as the 
Strength of the shock wave is increased. For air, the 
Mach Number limit is 1.89. By using gases with lower 
values of y such as carbon tetrachloride, this limit can 
be increased to about 4.5. It should be pointed out 
here that this Mach Number limitation does not apply 
in the steady-flow region behind the interface. How- 
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ever, in this region the test section temperature problem 
is of the same order as in a blow-down tunnel. 

The Mach Number Jimitation in the flow before the 
interface can be circumvented by employing a diverging 
nozzle before the test section. When the flow is super- 
sonic behind the shock wave, reflected waves are con- 
tinuously given off by the shock wave as it travels 
through the diverging section. These reflected waves 
coalesce to form a secondary shock wave. When the 
strength of the main shock is correctly chosen, this 
secondary wave is swept downstream, and, after its 
passage, steady flow is established in the test section. 
The Mach Number in the test section is then deter- 
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Fic. 1. Wave diagram of the transient following the passage 
of a shock wave through a diverging channel, computed by the 
method of characteristics. The air is initially at rest and the 
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initial Mach Number of the shock is 3. 
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TEST SECTION MACH NUMBER 

Fic. 2. Comparison between initial temperatures and test 
section temperatures for both blow-down tunnel and shock-tube 
tunnel. In the shock-tube calculations, the back pressure was 
assumed to be equal to the initial pressure in the low-pressure 
chamber. 
mined by the entrance Mach Number and the area 
ratio between the low-pressure chamber and the test 
section. Fig. 1 is a typical wave diagram, obtained by 
the method of characteristics,*® of the wave pattern 
accompanying a shock wave with steady supersonic 
flow behind it entering a diverging test section. 

It is an important property of this shock tube that the 
gas is heated by shock compression prior to the expan- 
sion, and this heating is more than sufficient to eliminate 
low test section temperature problem. This is brought 
out in Fig. 2, in which the ratio 7 of the static test sec- 
tion temperature to the temperature of the air at rest 
before acceleration to hypersonic flow is plotted for both 
this shock tube and the blow-down tunnel. 

In order to study further this method of flow genera- 
tion, a small shock tube was set up at C.A.L., and a 
series of exploratory tests were carried out using air in 
the low-pressure chamber and helium in the high-pres- 
sure chamber. The method of instrumentation em 
ployed is described in reference 10. Shock waves were 
fired into an expanding test section designed for an 


DECEMBER, 1951 
approximate Mach Number of 4.3. When the flow was 
supersonic behind the shock, the secondary wave was 
observed. This wave can be seen behind the main shock 
in Fig. 3, which shows the shock wave as it enters the 
test section. This picture also shows separation that 
occurred upstream of the secondary shock. This 
separation followed the secondary shock through the 
test section. Steady flow at Mach Number 4.2 was ob- 
served after the passage of this region of separation, 
Even the relatively low Mach Number of this explora- 
tory experiment was much higher than the Mach Num 
ber that can be created behind a shock front in air ina 
shock tube of constant area. Fig. 4 shows the bow 
wave from a wedge after the steady flow is established. 
Further studies into this method of flow generation are 
being carried out and will be published upon comple- 
tion. 
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PRIMARY SHOCK 


Fic. 3. Schlieren pictures of a shock wave during passage 
through a diverging channel. The initial strength of this wave 
was sufficient to produce supersonic flow at the entrance to the 
diverging section. Both the primary and secondary shock waves 
may be seen. 





Fic. 4. Steady flow established in test channel showing bow 
wave from test airfoil 
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On the Deflection of Swept Cantilevered 
Surfaces’ 


HAROLD C. MARTIN? ann HEMEN J. GURSAHANEYt 
Unwersity of Washington 


SUMMARY 


An experimental investigation is carried out for the purpose of 
determining the applicability of elementary theory in predicting 
the behavior of swept cantilevered surfaces. Attention is focused 
on swept plates, and primary conside1ation is given to deflections, 
although some information on stresses and dynamic character 
istics is also included. 

Plates with and without taper in plan form are tested at various 
angles of sweep. It is shown experimentally that, except in the 
neighborhood of the clamped edge, an axis exists which does not 
deflect under application of a pure twisting moment to the plate. 
Replacing the plate by an equivalent beam, located along this 
“apparent” elastic axis, permits conventional methods to be em- 
ployed for calculating deflections, bending frequencies, and so on. 
Such calculated values are shown to be in good agreement with 
experimental data, while differing considerably from results based 
on the conventional elastic axis. 


INTRODUCTION 


_ DESIGN OF SWEPT WINGS for high-speed aircraft 
and missiles will generally involve aeroelastic con- 
siderations. In other words, the wing must be treated 
as an elastic structure, and its deflection must be taken 
into account when computing aerodynamic character- 
istics. ' 

Various methods have been proposed for calculating 
the deflection of a swept wing under an arbitrary load- 
ing. The simplest of these is based on elementary 
beam theory and actually represents an application of 
methods developed for the unswept wing to the case of 
the swept wing. Essentially, such a procedure consists 
of postulating the existence of an elastic axis and then 
replacing the actual wing by an equivalent beam located 
along this axis. In the unswept wing, the location of 
the elastic axis depends only on the cross-sectional 
structural geometry. For either the unswept or swept 
wing, it is considered to be that axis that does not de- 
flect when a pure twisting moment is applied to the 
wing tip. (The twisting moment vector is here under- 
stood to agree in direction with the elastic axis. ) 

If an elastic axis does exist, bending and twisting may 


be considered separately. Since the problem may be 
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assumed linear, the deflection due to each separate load 
ing may then be added to give the final resultant de- 
flection. This procedure is relatively simple; there- 
fore, if also applicable for the swept wing, it would pro- 
vide a familiar process for dealing with the structural 
problems of such a wing. 

In this paper it is assumed that the essential influence 
of sweep on deflections will be of such a nature as not 
to be fundamentally altered if the actual wing struc- 
ture is replaced by a flat plate. The differences be- 
tween the actual wing and plate, particularly at the 
clamped end, cannot be stated with certainty. It can 
only be assumed that, in principle, the influence of 
sweep on deflections will be qualitatively the same in 
both cases. 

Actually, there are several reasons for first consider- 
One is that the mathematical 
Another 


ing the swept plate. 
formulation of the plate problem is known. 
is that the plate is simpler to obtain and test than an 
actual wing. 

A previous paper? by the authors has shown that for 
the swept rectangular plate an axis exists which cor- 
responds to the conventional elastic axis of the un- 
swept plate. The location of this axis depends on the 
angle of sweep, as well as on the cross-sectional proper- 
ties. Replacing the plate by an equivalent beam lo- 
cated along this apparent elastic axis will enable deflec- 
tions to be calculated with good accuracy by means of 
elementary theory. Similar calculations based on the 
so-called conventional elastic axis give, on the other 
hand, results that are often in error by 15 per cent or 
more. 

In the neighborhood of the clamped edge, conditions 
are complex and cannot be satisfactorily discussed on 
the basis of simple theory. Here, however, deflections 
are small, leading to negligible inaccuracies in the over- 
all aeroelastic calculations. 

This paper briefly discusses the experimental evi- 
dence leading to the concept of the apparent elastic axis 
for the swept plate of rectangular plan form. At the 
same time, a physical picture of the essential effect of 
Com- 
parison between experimental and calculated deflec- 
tions is presented for the swept rectangular plate sub- 
jected to a uniformly distributed load. Calculated 
values are based on both the conventional and apparent 


the plate root region on the deflection is given. 


elastic axes for further comparison. 
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With this preliminary information in mind, a detailed 
investigation is carried out for the swept plate of tri- 
angular plan form. It is shown that, for these plates 
also, an apparent elastic exists which permits the plate 
to be replaced by an equivalent beam. Experimental 
data are compared with calculated results for different 
loadings on plates having several angles of sweep. 
Finally, the first three natural bending modes of vibra- 
tion are investigated for the swept rectangular plate. 
It is shown that both the natural frequencies and mode 
shapes may be closely predicted by again considering 
the plate as a beam located at the apparent elastic axis. 


THE SWEPT RECTANGULAR PLATE 


In reference 2 it is shown that the deflection of a swept 
rectangular plate may be calculated by elementary 
theory, provided the plate is replaced by an equivalent 
beam located at the apparent elastic axis. This axis is 
parallel to the spanwise centerline* of the plate but is 
shifted toward the trailing edge as shown in Fig. 1. 
Its location cannot be determined from simple theory 
and would be difficult to calculate from basic plate 
theory. Fig. 1 therefore represents experimental data 
(taken from reference 2). 


* The elastic axis of the unswept plate—referred to in this 
paper as the conventional elastic axis. 







































APPARENT Q 
Ak ELASTIC Axis i 
‘ = 
e T T 
o EXPERIMENTAL POINTS 
/ 
ft 
8 
7 y 
6 
Q% = 0725-0225 Cos 2V 
5 
fe) 20 40 60 


Y- DEGREES 
Location of the apparent elastic axis (swept rectangular 
plate). 


Fic. 1. 


AERONAUTICAL 





SCIENCES DECEMBER, 1951 

When a pure twisting moment is applied at the tip of 
a swept rectangular plate, deflection measurements 
show that the conventional elastic axis deflects. By 
definition, the apparent elastic axis does not deflect 
under such a moment. It is also significant that the 
deflection of the conventional elastic axis is due to in- 
duced bending in the root portion of the plate—that is, 
the pure twisting moment applied at the tip induces bend- 
ing near the clamped edge. As a result, the deflection 
curve for the conventional elastic axis shows curvature 
near the clamped edge but a constant slope in the outer 
span region. 

These experimental observations are further de- 
veloped in reference 2. Various loadings, in addition to 
that mentioned above, are applied to the plate in order 
to verify calculations based on the apparent elastic axis. 
Of these, the uniform loading is of greatest interest to 
the airplane designer, since it imparts variable bend- 
ing and twisting moments to the plate in a manner 
analogous to that caused by flight loads. Further- 
more, the uniform load has another interesting feature 
in this case—namely, its center of pressure coincides 
with the conventional elastic axis. Therefore, deflection 
calculations based on this axis will include bending 
only, while those based on the apparent elastic axis will 
include twisting, as well as bending. 

Experimental and calculated data may be conven- 
iently compared by forming the ratio 


— experimental result (1) 
~ calculated result 
Deviations of k from unity then measure the error in 
calculated values. 

Table 1 gives the k-factors for a rectangular plate of 
aspect ratio 7 and having an angle of sweep, Y = 40 
The total load on the plate (uniformly distributed) is 
43.9 Ibs. Calculated values are given for both the con- 
ventional and apparent elastic axes. In either case, the 
uniform beam loading (pounds per inch) is the beam 
length divided into the total loading on the plate. 

Stresses were not experimentally determined during 
the course of this investigation. Such data are, how- 
ever, presented in reference 3 for swept rectangular 
plates. Calculations given in reference 2 show that, 
except near the clamped edge, bending and shear stresses 
may be calculated with reasonable accuracy by elemen- 
tary methods. Near the clamped edge discrepancies 
between theory and experiment can be expected to be 
large, since simple theory cannot adequately describe 
the complex situation existing in that portion of the 
plate. 


THE SWEPT TRIANGULAR PLATE 


Introduction 


Actual cantilever lifting surfaces will generally be 
tapered, and therefore the effects of taper on deflections 
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TABLE 1 


k (Based on Apparent 
- - Elastic Axis) 


Leading Trailing 


Station Leading Spanwise Spanwise 
&/l’ edge edge edge edge 
o.3 2.150 1.939 1.1hi 
0.2 ; a : 1.135 1.073 1.029 
0.3 0.837 0.830 0.815 0.970 0.982 0.984 
0.5 : ; 0.949 0.955 0.950 
0.7 0.895 0.860 0.856 0.946 0.943 0.945 
0.9 0.982 0.971 0.966 
1.0 0.916 0.905 0.897 0.9838 0.971 0.964 


are studied next. The primary purpose of such an in- 
vestigation will be to see whether or not a straight- 
line apparent elastic axis exists for the tapered canti- 
lever. 

Since the triangular plate (uniform thickness) repre- 
sents the limiting case of complete taper (in the sense 
that the rectangular plate represents complete ab- 
sence of taper), it appears logical to select triangular 
plates for investigation. As will be seen later, results 
obtained from plates of triangular plan form appar- 
ently have significance for the arbitrarily tapered cases 
also. 

It is not obvious that an apparent elastic axis, which 
is also a straight line, exists for the swept triangular 
plate. For example, such an axis implies that a twist- 
ing moment vector, applied so that its direction coincides 
with the axis, would leave the plate tip undeflected. 
This is so because the apparent elastic axis must include 
the plate tip and by definition must remain undeflected 
under application of a twisting moment oriented as 
above. (It is to be understood that these remarks apply 
essentially to the outer span portion of plates having a 
moderately high aspect ratio.) 

Assuming that the straight-line apparent elastic does 
exist for the swept triangular plate, it is clear that its 
location must be purely a function of the structural 
geometry of the plate. In other words, this axis is 
unique and its location depends on the physical struc- 
ture only, being independent of the applied loading. 

The above ideas are useful in experimentally locating 


the desired apparent elastic axis. 


Test Specimens and Loading 

Typical test plates are illustrated in Fig. 2. The 
angle of sweep, y, is measured to the spanwise plate 
centerline, 0'’A. This line becomes the conventional 
elastic axis when y = 0. 

In order to be sure that a considerable portion of the 
plate area remains remote from the influence of the 
clamped edge, test specimens were given an aspect 
ratio of 10. The definition of aspect ratio (A.R.) 
used here is similar to that for the airplane or 


A.R. = (span)?/area = (2b’)?/b’c’ = 4b'/c’ 


Descriptive data for two test plates are given in 


Table 2. 


Actual testing was performed by clamping the plates 
between rigid structural members. Duplicability of 
readings for “‘load-off’’ and ‘‘load-on”’ conditions indi- 
cated absence of slip. Deflections were measured with 
a standard dial gage supported on a smooth surface 
beneath the test specimen; this surface was supported 
independent of the specimen and the structure holding 
the specimen. 

Loadings used were dead weights applied as simply as 
possible. For example, suppose a twisting moment of 
vector direction coinciding with 0’A in Fig. 2 were to be 
applied near the plate tip. Such a loading is simply 
obtained by mounting a light bracket onto the plate 
near the tip and normal to 0’A. A given dead weight 
applied at two different points on this bracket then pro- 
duces a twisting moment about 0’A of magnitude equal 
to the weight times the distance between its two points 
of application. The difference in deflection readings, 
taken for each position of the load, then represents the 
deflection due to the twisting moment only. Similar 
means can be used to apply a pure bending moment to 


the plate. 


Pure Twist Loading 


A pure twisting moment is defined with the apparent 
elastic axis in mind. Assuming again such am axis to be 
an invariant for a given plate at a given angle of sweep, 
this axis will, by definition, remain undeflected under 
application of a twisting moment of vector direction 
parallel to its length. 








b= % 














Fic. 2. Swept triangular plate. 
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TABLE 2 ; 
v b’ c’ h E G 7 
Specimen ( Deg.) (In.) (In.) A.R. (In.) (Lbs. per Sq.In.) (Lbs. per Sq.In.) 
A 25 30 12 10 5/16 10.3 * 108 3.86 X 108 
B 43.8 25 10 10 5/16 10.3 X 108 3.86 X 108 


TABLE 3 
y = 25°, M, = 162.5 In.Lbs. (Pure Twisting Moment Applied at £// = 0.85) 


Calculated 
—- Deflections 


E 0 Leading Trailing 
E/1 (In.) (rad.) edge edge 
0.1 3.23 0.001331 —0.0078 0.0044 
0.2 6.46 0.002817 —0.0147 0.0817 
0.3 9.69 0.00453 —0.0197 0.0118 
0.5 16.15 0.00876 —0.0279 0.0152 
0.7 22.61 0.01522 —0.0283 0.0146 
0.8 25.84 0.02038 —0.0236 0.0120 


A logical first step in locating this axis is to so orient 
the twist vector that the plate tip remains undeflected. 
When this has been achieved, one or two other spanwise 
points of zero deflection may likewise be found. A 
straight line from the tip through these other points 
of zero deflection tentatively locates the apparent elas- 
tic axis. 

This location of the apparent elastic axis should then 
be carefully checked. Calling the spanwise coordinate 
along the tentative axis direction, £, the leading- and 
trailing-edge deflections may be plotted as in Fig. 3. 
These deflections are for the twisting moment that gave 
zero tip deflection. The dotted portions indicate plate 
area blocked out by the load support bracket. 

From Fig. 3 it is a straightforward matter to plot 
chordwise views of the deflection at selected stations 
along the span. Such sections will be normal to the 
apparent elastic axis if the leading-edge and trailing- 
edge deflections are taken from Fig. 3 at a given sta- 
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Fic. 3. Deflection due to pure twisting moment (y = 25°, AZ; = 


162.5 in.Ibs. ). 


Experimental 


Deflections - — k [See Eq. (1)] 





Leading Trailing Leading Trailing 
edge edge edge edge 
—0.0064 0.0019 0.821 0.482 
—(0.0125 0.0067 0.850 0.850 
—0.0185 0.0118 0.939 0.958 
— 0.0266 0.0153 0.954 1.006 
—0.0276 0.0151 0.989 1.027 
—().0240 0.0123 1.016 1.025 


Although not shown, such plots were made 
For each 


tion, &/1. 
for sections at £/] = 0.2, 0.3, 0.5, and 0.7. 
such plot, the point of zero deflection was noted and 
plotted on a plan view of the plate, as in Fig. 4. It is 
readily seen that the curve through these points of zero 
deflection and including the plate tip is straight line. 
This straight line is continued into the plate clamped 
Line OA, Fig. 4, is then the apparent 
In carrying out the above procedure, the 


edge at 0. 
elastic axis. 
orientation of the twisting moment vector may have to 
be altered slightly in order to obtain the straight line, 
0A, shown in Fig. 4. 

The process just carried out may now be repeated 
for plates of the same A.R. but having different angles 
of sweep. Two such plates give data as shown in Fig. 
5; also, in this figure, is the result of the apparent elastic 
axis determination for a plate of A.R. = 7 and sweep, 
y = 40°. The results for all three plates indicate that 
the intersection of the straight-line apparent elastic axis 
with the plate clamped edged may be represented by 
the linear relationship 

a/c’ = 0.5 + (w°/150) (2) 
where a and c’ are shown in Fig. 5. This figure also 
shows a plot of Eq. (2) and how it fits the experimental 
data. 

Having established the location of the apparent elas 
tic axis for swept triangular plates, it now becomes a 
question of some interest as to how accurately deflec- 
tions may be calculated when the plate is reduced to a 
beam located along this axis. For the plate in pure 
twisting, the basic elementary equation for angle of 


twist is 
dé/dé = M,/GJ(é) (3) 
where 
6 = angle of twist of section normal to the ap- 
parent elastic axis (0) 
M, = applied twisting moment, in.Ibs. 
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DEFLECTION OF SWEPT 

J(é) = Jo(l — &)/land/ = length of apparent elastic 
axis, while J) is the section stiffness in 
torsion’ at the root (Coh*/3 for a rectangle 
with Cy/h > 10) 


Eq. (3) may be integrated at once to give 
6 = (M,l/GJo) log 1/(l — &) (4) 


By means of Eq. (4) the twist of any section normal to 
axis 0& may be calculated. Knowing 6 and making use 
of the fact that axis 0 does not deflect, the deflection, 
w, at any chordwise position may then easily be deter- 
This has been done for test plate A, and re- 
Again, experimental and 


mined. 
sults are given in Table 3. 
calculated deflections are compared on the basis of a 
k-factor found from Eq. (1). 

Agreement between theoretical and experimental re- 
sults is seen to be encouraging. 


Concentrated Load at Tip 

Since the apparent elastic axis passes through the 
plate tip, a concentrated load placed at this point should 
produce bending (only) of the equivalent beam. Such 
a loading, therefore, is useful in further checking the 
correctness and significance of the apparent elastic axis. 

In performing such tests, plates A and B of Table 2 
were subjected to concentrated loads placed at the plate 
tip of 13 and 8 lbs., respectively. The leading- and trail- 
ing-edge spanwise deflections for specimen A are shown 
Constant £// represents a specific section 
This illustration 


in Fig. 6. 
normal to the apparent elastic axis. 
clearly indicates that such chordwise sections deflect 
without twisting under the applied tip load. Near the 
clamped edge a small amount of induced twisting is seen 
to occur. A similar plot for specimen B duplicates the 
behavior already shown for specimen A. 

Once again, comparison between calculated and ex- 
perimental deflections is significant. The theoretical 
values are obtained by substituting a cantilever beam, 
located at the apparent elastic axis, for the triangular 
plate. The deflection at any point, &, is then given by 
elementary theory as 

w= Ple/2EIy (5) 
where /) is the moment of inertia of the root section 
and equals Cjh*/12. Table 4 gives experimental de- 
flections for specimen A and k-factors for both test 
plates A and B. 

Simple theory again proves to be accurate for calcu- 
lating deflections in that region of the plate removed 
from the immediate vicinity of the clamped edge. 


Uniformly Distributed Load 


As in the case of plates of rectangular plan form, the 
combined twisting and bending of the triangular plate 
subjected to actual flight loading is conveniently dupli- 
cated in its essential details by a uniformly distributed 


load. Furthermore, such a loading can be easily ap- 
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plied in the laboratory and handled without excessive 
labor in theoretical calculations. Consequently, plates 
of triangular planform are investigated for such a load- 
ing. 

The basic differential equation of beam bending 


d?w,/dé? = M,(é)/EI(é) (6) 


is used in calculating bending deflections. Twisting 
deflections may be calculated by applying Eq. (4) 
and by noting that /, is now a function of &. 

If g is the uniform pressure (pounds per square inch) 


on the plate, the bending moment will be given by 





















































810 JOURNAL OF THE AERONAUTICAL 
/ 
‘ 
07 p 
0.6 
o EXPERIMENTAL POINTS 

0.5 

= 
' 
z O4 
© 
Ee 
UU 
4 
u 03 " 
lJ 
a 
Oo 
< 02 
s 
= 
re) WA 
Oo al 
LEADING 
EOGE 
TRAILING EDGE 
Oo 0.2 0.4 06 08 1.0 
F/e 

Fic. 6. Deflection due to tip force (Y = 25°, P = 18 Ibs.). 


M,(é) = qC(l — £&)°/6 


where C is the chord length normal to the apparent 
elastic axis 0 at — and may be written as C,(/ — &)/l. 
Integrating Eq. (6) twice and imposing the boundary 
conditions that the deflection and slope of the equivalent 
beam vanish at — = 0 gives, for the bending deflection, 


wy = (Cog/72EIy) (2 — &)4 + 40 — I] (7) 


The total deflection is the sum of the deflections due 
to bending and twisting, here calculated separately. 
This is a matter of detail, since the necessary Eqs. (4) 
Results of such calculations 
The loading used 


and (7) are now available. 

are given in Table 5 for test plate A. 

is g = 0.1744 Ibs. per sq.in. 
Experimental values for the deflection are determined 


for this same loading. Actually, dry sand, contained 


SCIENCES—DECEMBER, 1951 
trailing edges of the plate, was used in loading the plate. 
It is interesting to note that any stiffening introduced 
by the cardboard was released by placing the dial gage 
at the plate tip and then cutting horizontal and vertical 
slits in the cardboard until deflections became unaffected 
by such action. The total change in deflection observed 
during this process was of the order of 5 per cent. 

Table 5 also lists the measured deflections and the cor- 
responding & factors comparing calculated and experi- 
mental results. The accuracy of the elementary theory 
is seen to be remarkably good. 

The results shown in Table 5 substantiate previous 
results in that, at least as far as deflection calculations 
are concerned, the swept triangular plate may be treated 
as a cantilever beam if located at the plate apparent 
elastic axis. 

Since the conventional and apparent elastic axes inter- 
sect at the plate tip for swept triangular plates, they 
are not separated in the outer span region to the same 
extent as in the case of rectangular plates. Conse- 
quently, the distance between the center of pressure of 
the uniformly distributed load and the apparent elastic 
axis is relatively less significant for the triangular plates. 
In other words, calculated results based on the conven- 
tional elastic axis should not deviate so greatly from 
those based on the apparent elastic axis as was pre- 
viously observed for the rectangular plates (see Table 
1). Table 6 indicates this to be true. Comparisons 
shown are for tip deflections only. 


Arbitrary Taper 

The arbitrary tapered plan form is of interest, since 
most swept lifting surfaces are of such design. For the 
plates discussed here, arbitrary taper will have the rec- 
tangular and triangular plates as limiting cases. It 
seems reasonable to conclude, therefore, that the ap- 
parent elastic axis will again be a straight line. 

As indicated by Eq. (2) and Fig. 5, the location of the 
apparent elastic axis for the triangular plate depends on 
sweep only—that is, is independent of semispan, 0’. 
This invariance is analogous to that previously ob- 
served for the rectangular plates. In this latter case, 
the location of the apparent elastic axis depended also 
on y only and was not affected by the plate semispan, 


between cardboard sides attached to the leading and 6’. Again, these remarks are understood to apply only 
TABLE 4 
(Concentrated Load at £//] = 0.994 for Each Specimen ) 
———_—_—_—____——_ Plate A —_— - Plate B —- 
w (experimental) w (experimental) ———k — k 
w (calculated ) leading trailing Leading Trailing Leading Trailing 
E/1 (In.) edge edge edge edge edge edge 
0.1 0.0078 0.0145 0.0090 1.858 1.154 1.191 0.888 
0.2 0.0312 0.0400 0.0360 1.278 1.154 0.991 0.959 
0.3 0.0702 0.0800 0.0775 1.140 1.104 0.954 0.954 
0.5 0.1951 0.2050 0.2050 1.051 1.051 0.957 0.957 
0.7 0.3824 0.3880 0.3880 1.015 1.015 0.958 0.958 
0.9 0.€104 0.6200 0.6200 1.015 1.015 0.963 0.963 
0.969 (A) 0.7328 0.7110 0.7110 0.970 0.970 
0.945 (B) 0.974 0.974 
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TABLE 5 
gq = 0.1744 Lbs. per Sq. In, y = 25°, 1] = 32.3 In., Jo = 0.02708 In.4 
7 Total Calculated Experimental 
———Deflection ~ - Deflection — — k ~ 
We @ xX 108 Leading Trailing Leading Trailing Leading Trailing 
t/l (In.) (Rad.) edge edge edge edge edge edge 
0.1 0.0012 0.179 ().0022 0.0001 0.0078 0.0060 
0.2 0.0211 0.249 0.0224 0.0204 0.0275 0.0230 1.228 1.128 
0.3 0.0443 0.269 0.0455 0.0436 0.0500 0.0465 1.099 1.066 
0.5 0.1070 0.186 0.1076 0. 1067 0.1130 0.1130 1.051 1.059 
0.7 0.1820 0.076 0.1821 0.1819 0). 1860 0.1860 1.020 1.022 
0.9 0.2618 0.054 (0). 2618 ().2618 0.2562 (0). 2562 0.970 0.975 
1.0 0.3020 0 0.3020 0.3020 0.2920 0.2920 0.966 0.966 
to plates of sufficiently high A.R. to put a reasonable where 
proportion of the plate area beyond the immediate 
, : a eer ee f. = cycles per sec. 
neighborhood of the clamped edge. The lower limit of wd 4 : ; . 
: . EI = flexural stiffness in bending 
such aspect ratios cannot be stated from data on hand. ating: , 
l =} length of beam, which is the length of the ap- 
A first approximation for the apparent elastic axis parent elastic axis for the plate (31.5 in.) 
location of the arbitrarily tapered plate might, perhaps, “ = mass/unit length of plate = 322.6 * 10~® Ibs. 
be obtained by extending the plate to a triangle and sec.2/in.? 
then using Eq. (2). This has not been checked experi @, = numerical coefficient dependent on mode. 
mentally, however. n = 1, 2, 3, represent the first three bend- 
ing modes and for these a, = 3.52, a = 
DYNAMIC CHARACTERISTICS 22.4, ds = 61.7 
In the above, » = 1 represents the fundamental mode 


Introduction 

A preliminary investigation indicated that at least 
some of the basic dynamic characteristics of the swept 
plate could be explained on the basis of elementary 
theory. Experimental data for the free vibrations in 
bending of the swept rectangular plate were compared 
with calculated values. Quantities measured were: 
first three natural bending frequencies, the shape of the 
fundamental bending mode, and the nodal lines on the 
plate for the second and third bending modes. 


Experimental Method 

Tests were conducted on a rectangular swept plate 
having y = 50° and A.R. = 7. Excitation of the nat- 
ural frequencies in bending was achieved by means of 
an electromagnetic shaker* located on the apparent 
elastic axis at the plate tip. Deflections during vibra- 
tion were measured with electrical pickup units de- 
veloped at Boeing. 

Nodal lines occurring at the second and third natural 
modes were obtained by sprinkling sand on the upper 
plate surface. This sand moved to clearly defined 
lines as the natural modes were excited. 


Natural Frequencies in Bending 

Elementary beam theory gives the following expres- 
sion’ for the natural bending frequency of the canti- 
lever: 


f, = (d,/2e) WEI /pl' (8) 


* Vibration equipment was generously furnished by the Boeing 
Airplane Company. 


in bending. 

In Table 7 calculated natural frequencies from Eq. 
(S) are compared with experimentally determined 
values. It is seen from this table that results obtained 
by using the conventional elastic axis are in consider- 
able error compared to those obtained by use of the ap- 
parent elastic axis. 


Nodal Lines 


The nodal lines as obtained by sprinkling sand on the 
plate are shown by the dotted lines in Fig. 7. These 
lines were found to be inclined to the apparent elastic 


axis as shown. This inclination is probably due to the 


TABLE 6 


k (Concentrated k (Uniformly 
-——Tip Load)——. —Distributed Load)—~ 


Con- 
Appar- ven- Appar- Con- 
. ent tional ent ven- 
Speci- y elastic elastic elastic tional 
men (Deg.) axis axis axis axis 
A 25 0.958 0.909 0.966 0.920 
B 43.8 0.962 0.883 
TABLE 7 
——Cycles per Min. 
Calculation k 
using Using Using 
Natural apparent apparent conventional 
Bending Experi- elastic elastic elastic 
Mode mental axis axis axis 
1 624 604 1.032 0.910 
2 3,810 3,840 0.992 0.872 
3 10,420 10,500 0.994 0.874 
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Fic. 7. Nodal lines in natural bending modes of vibration 
(y = 50°, A.R. = 7). 
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A.R. = 7). 


complex induced twisting taking place in the root region 
of the plate. 

For the cantilever beam in bending, the nodal points 
are well known and have been tabulated for easy refer- 
These points representing elementary theory are 
It is clearly seen that the nodal 


ence.°® 
also shown on Fig. 7. 
lines for the plate cross the apparent elastic axis at 
virtually the points given by beam theory. 
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TABLE 8 a 

— —w/Wip . : 
E/1 Experimental Calculated k 
0.15 ine 0.0274 as 
0.30 0.1380 0.1089 1.16 
0.45 0.2332 0.2394 0.98 
0.60 0.4230 0.4118 1.03 
0.80 0.6960 0.6904 1.01 
0.95 0.9322 0.9342 1.00 
1.00 1.0000 1.0000 1.00 


First (or Fundamental) Bending Mode 


The deflected shape of the cantilever beam vibrating 
in its fundamental mode is given by elementary theory 
as? 


W/Wtip = 1 — cos (ré/2/) (9) 


where w is the deflection and £ is the spanwise coordi- 
nate measured from the clamped end. The deflection 
ratios as calculated from Eq. (9) for a number of span- 
wise locations, £//, are given in Table 8. These are 
based on the apparent elastic axis—that is, / = 31.5 
for the plate having y = 50° and A.R. = 7. 
Experimental deflections at these same spanwise 
points on the apparent elastic axis were also deter- 
mined. 
between theory and experiment is seen to be excellent. 


These are also given in Table 8. Agreement 
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Dynamic Effects in Rotor Blade Bending’ 


A. H. FLAX ann L. GOLAND# 


Cornell Aeronautical Laboratory, Inc. 


ABSTRACT 


The problem of helicopter rotor-blade bending moments for the 
case in which the blade is near resonance with one of the periodic 
aerodynamic forcing harmonics is considered. The additional 
inertia and aerodynamic damping forces due to blade bending 
deflection are included in the genera! equations for blade bend- 
ing. Two methods of analysis are presented: The first is an ap- 
proximate method based on the use of an amplification factor on 
the static bending moments of the blade; the second is a tabular 
solution of the differential equation. The tabular procedure is 
rather laborious and time-consuming, while the amplification 
factor method is relatively simple. Comparison of results ob- 
tained by the two methods for a typical case indicates fairly good 


agreement. 


INTRODUCTION 


I A PREVIOUS PAPER by one of the authors,! the vari- 
ous available methods for the solution of the prob- 
lem of rotor-blade bending were reviewed and com- 
pared, and certain extensions of these methods were 
presented. All of this work involved the assumption 
that the additional inertia and aerodynamic forces due 
to blade deflection could be neglected. It was pointed 
out that this presupposed that the periodic aerody- 
namic forces were not near resonance with one of the 
natural bending frequencies of the blade. The analysis 
of the C30 autogyro rotor blade by Hufton, Nutt, Bigg, 
and Beavan? was cited as evidence that such an ap- 
proach might be valid for typical rotor blades. More 
recently, however, Johnson and Mayne* have shown 
that, for a particular but typical helicopter rotor blade, 
there was in fact a resonance of the second harmonic 
aerodynamic forcing function with the fundamental bend- 
ing frequency of the rotating blade. This condition 
obviously necessitates a consideration of the additional 
blade forces caused by blade deflection. Johnson and 
Mayne, in their work, took into consideration the ad- 
ditional inertia forces due to blade deflection but did not 
attempt to account for the additional aerodynamic 
forces. They solved the problem by means of a lengthy 
tabular method and then proposed to treat this problem 
practically by multiplying the second harmonic ‘‘static’’ 
blade moments by the amplification factor for a simple 
undamped single degree of freedom system. Hohen- 

Presented at the Rotating Wing Aircraft Session, Nineteenth 
Annual Meeting, I.A.S., New York, January 29-February 1, 
1951. 

* This investigation was carried out under sponsorship of the 
Rotary Wing Branch, Air Materiel Command, to which the 
authors are appreciative for permission to publish the results. 

7 Head, Aerodynamic Research Department. 

t Associate Research Engineer, Aero-Mechanics Department. 


emser‘ had previously indicated the possibility of deal- 
ing with the problem in this way but also pointed out 
that the additional aerodynamic forces that are prin- 
cipally damping forces might alter the results con- 
siderably. It was also observed in certain applications 
of the Johnson and Mayne amplification factor to 
cases where the forcing harmonics were extremely close 
to coincidence with the blade natural frequency that 
unreasonably large stresses were computed. Obvi- 
ously, the aerodynamic damping would limit the ampli- 
fication factors in these cases. The present study was 
undertaken to evaluate the effects of such damping on 
the blade bending moments. 

Two methods of analysis are developed—an approxi- 
mate amplification factor method and a more precise 
tabular method. The approximate method presented 
herein is essentially that presented in reference 5, an 
earlier report on this work. The tabular method is 
rather difficult to apply, since it requires lengthy and 
accurate computations; the approximate analysis offers 
a rational and feasible method by which solutions are 
easily obtained without any tedious calculations. In 
fact, application of the approximate method involves 
only a small amount of additional work over that re- 
quired for the Johnson and Mayne amplification factor 
method. Outlines of procedures for the application of 
the two methods, which may be used independently of 
the theoretical development, are presented. 

As an application of the methods, a typical 19-ft. 
rotor blade is analyzed for a condition of forward level 
flight (120 m.p.h., 252 r.p.m.). 
sented between the results obtained by the two methods 
developed in this paper-and by the two methods sug- 
It is found that inclusion of the 


Comparisons are pre- 


gested in reference 3. 
aerodynamic damping considerably reduces the ampli- 
fication of the blade moments and that, in addition, it 
introduces a phase lag between the blade position at 
maximum loading and maximum moment. 

Agreement between the magnitudes of the second 
harmonic bending moments obtained by the approxi- 
mate and more exact methods is very good. The phase 
lags do not agree so well, but the approximate values are 
satisfactory for most purposes. After this work had 
been completed, the authors’ attention was called to a 
report by de Guillenchmidt,® which has been translated 
from the French by the N.A.C.A., presenting a method 
similar in many respects to the approximate method. 
There are, however, some significant differences be- 
tween the two methods, and it is believed that the way 
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in which damping is introduced in this paper is more 
accurate. Further, the approximate method presented 
herein makes use of the more accurate tabular static 
solution, while reference 6 is based entirely on the ap- 


proximate analysis. 


NOTATION 


distance from the center of rotation to a blade sec- 


r = 
tion 

(r) = symbol denoting that a quantity is a function of r 

aor) = aerodynamic slope of the lift curve of the blade air- 
foil section 

I(r) = sectional moment of inertia 

E(r) = modulus of elasticity of blade material 

c(r) = blade chord 

m(r) = mass of blade per unit length 

R = length of the blade 

X = length of interval between blade stations 

n = number of intervals into which the blade is divided = 
R/X 

a(r = aerodynamic damping coefficient 

Z(r deflection of the elastic blade from the corresponding 
position of the blade considered as rigid and 
hinged at the rotor hub 

F(r) = total centrifugal force at any station r of the blade, 
obtained by summing up all incremental centrif- 
ugal forces between the point 7 and the tip 

p(r) = axial loading (centrifugal force) per unit length of the 
blade 

S(r) = the net shear acting on a rigid blade hinged at the 
rotor hub, obtained by summing up all incremental 
loads between the point 7 and the tip (positive 
downward ) 

D(r) = the aerodynamic damping shear due to the elastic 
blade deflection, obtained by summing up all in- 
cremental aerodynamic damping forces between 
the point 7 and the tip (positive downward ) 

T(r) = the “inertia force’’ shear due to the elastic blade de- 
flection, obtained by the summation of the incre- 
mental inertia forces between the point 7 and the 
tip (positive downward ) 

H(r) = total shear at any station 7 of the elastic blade (posi- 
tive downward ) 

V = helicopter forward velocity 

m = advanceratio = V/QR 

M(r) = bending moment at any station 7, represented in the 
amplification factor method by 
M = M + Miu sin Qt + Min cos Qf + 

Myy sin 20t + My cos 20t +... 
and in the tabular method by 
M = M' cos kQt + M" sin kQt where k = 
ay ee: 

M'(r) = “‘static’? bending moment obtained by neglecting 
the damping and inertia force loadings due to 
blade bending: 

My + My’ sin Qt + Mint’ cos Qt + Miy’ sin 20t + 
My’ cos 2Qt +... 

th = transverse aerodynamic and inertia loading per unit 
length acting on the blade (positive upward ) 

be = transverse aerodynamic loading per unit length on a 


rigid blade hinged at the rotor hub 
/ transverse loading due to aerodynamic, inertia, and 
centrifugal forces on a rigid blade hinged at the 


rotor hub 
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tk = kth harmonic transverse loading per unit length on 
the rigid blade hinged at the rotor hub 

Pn = dimensionless deflection of rotating blade in mth nor- 
mal bending deflection mode 

Wn, = frequency of mth normal bending deflection mode 

A = generalized deflection coordinate of nth normal mode 

y = azimuth position of blade measured from downwind 
position 

p = mass density of air 

B = hinged rigid blade flapping angle = a) — a cos y — 
b, sin Y — azcos 2y — besin 2y, ete. 

ak = cosine component of kth harmonic flapping angle 
(see B) 

b; = sine component of kth harmonic flapping angle 
(see B) 

t = time 

k = order of harmonic considered 

Q = rotational velocity of the rotor blade 

i = letter designating a particular blade station 

i = “running” index used in summations 

THE AMPLIFICATION FactoR METHOD INCLUDING 


AERODYNAMIC DAMPING 


Theoretical Development 


The basic equation for blade bending in the plane of 
flapping has been developed! in the form 


oO” (21 *) Flr) 0°Z 4 OZ ; ' 
or? \ Or? ” or? co ot PB) 


where 


4i=t,- ( f =) At ) 
mir + + 2 
L a #2 2 (<) 


Here, Af is the increment of aerodynamic loading 
caused by the deflection of the blade. This term may 
easily be derived by consideration of the increment of 
transverse velocity of the air relative to the blade, AV’z, 


produced by the bending motion of the blade. This 
leads to 
p ace AV; 
At = —d = c(Qr + V sin yp)? See 
2 (Qr + V sin y) 
Now, since 
AVz = 0Z/ot 
Al ‘ o+ Fae = 
= —M —Cc | (Mir sin y) 
> ar | 


and, consequently, At may be written in the form 


At = —a(r, ¥) (0Z/Ot) 


The use of the absolute value of the velocity at a sec- 
tion of blade rather than the algebraic value of the 
velocity is made necessary by the fact that Af remains a 
damping load even in the region of reversed flow on the 
retreating blade. It will be noted that this damping 
loading has been calculated on the basis of quasi- 
stationary aerodynamic theory. This is in line with the 
usual assumption made in calculating loads of rigid 
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DYNAMIC EFFECTS OF 
blades, and it appears to be satisfactory for the lower 
loading harmonics on conventional designs. To illus- 
trate, the reduced frequency, wc/2V, of a typical blade 
(R = 19 ft., c = 1 ft.) at its */, radius in second har- 
monic oscillation is about 0.07. According to the two- 
dimensional theory of nonsteady airfoil motion,’ the 
peak sinusoidal lift is 95 per cent of its quasi-stationary 
value for an airfoil oscillating at this value of the re- 
duced frequency, but it has a phase lag of about 10° of 
the vibration cycle or, in this case, 5° of the rotor azi- 
muth angle. Biot and Boehnlein’s studies of finite 
aspect ratio effects on unsteady airfoil theory show a 
considerable reduction of this small phase lag due to 
three-dimensional aerodynamic effects. While these 
calculations are not strictly applicable to rotor blades, 
it appears that the orders of magnitude indicated 
should be correct, thus justifying the neglect of un- 
steady aerodynamic effects, at least for the lower fore- 
ing harmonics. It should be noted that Isaacs” ex- 
tension of the two-dimensional theory to a flow of 
varying stream velocity, simulating the condition of a 
rotor blade, shows a similar small order of magnitude for 
unsteady effects. These arguments concerning the 
order of magnitude of the two-dimensional unsteady 
flow effects on rotor-blade forces, of course, become less 
and less valid as higher loading harmonics, with cor- 
respondingly higher frequencies, are considered. 


According to the usual elementary blade element 
theory of rotor loading, the contribution of the mth 
harmonic is proportional to wu". Thus, it might be ex- 
pected that the high harmonics would not make an 
important contribution to the blade bending stresses. 
Recent experimental work carried out at Cornell Aero- 
nautical Laboratory, Inc., under Air Force contract 
has indicated,“ however, that this may not always be 
true and that proper attention must be given to the con- 
tribution of the higher harmonics, although the loading 
is apparently not predictable from the usual blade ele- 
ment theory with simple downwash assumptions. 


Another point concerning unsteady flow effects must 
also be considered. At fairly high forward speeds, 
the wake from the rotor trails backward along the line 
of flight in much the same manner as a wing wake, so 
that comparative arguments may be considered to have 
some validity. At low forward speeds, however, 
where the flight velocities and downwash velocities 
through the rotor are of the same order of magnitude, 
the wake moves downward and backward away from 
the rotor at velocities that are small compared to the 
blade velocity, so that the geometrical configuration 
of the wake, which is the cause of the unsteady flow 
effects, is radically different from that considered in 
either two-dimensional or three-dimensional airfoil 
Thus, all of the above arguments must be 


theory. 
These 


considered invalid for low flight velocities. 
considerations concerning the aerodynamic forces on 
the rotor should be taken carefully into account in de- 
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ciding upon the degree of complexity to be tolerated in 
a rotor-blade analysis based on such approximate 
aerodynamic loadings. 

A further complication arises from the possible effect 
of torsional deflection of the blade. If this is appre- 
ciable and if the blade is not suitably mass balanced, 
the blade may approach a flutter condition. If this 
occurs, the aerodynamic damping will be considerably 
reduced, and the analyses given in this paper will be 
invalid. For purposes of this paper, it is assumed that 
the blades are either extremely stiff torsionally or well 
balanced and far removed from any flutter condition. 

It is seen from Eq. (3) that the damping load, Af, varies 

Inclusion of this effect in the 
to considerable complication, 


with azimuth position. 
analysis would lead 
since Eq. (1) would then have coefficients varying sinu- 
soidally with time. The exact solution would then re- 
quire Mathieu functions, and even an approximate 
solution would involve a considerable increase in the 
labor required. For this reason, the present analysis is 
made assuming a constant with y. Bennett," in the 
study of the stability of flapping hinged rigid blades, 
encountered a similar situation. His results indicate 
that the effect of this variation in the coefficients of the 
differential equation of motion is not sufficient to cause 
any blade instabilities even in the extreme case of an 
advance ratio, wu, of unity. For u < 0.5, which is typi- 
cal of current helicopters, the effect is even less severe, 
giving some assurance that fundamental changes in 
the character of the blade motion do not occur. A con- 
servative approach would be to use the minimum value 
of |Qr + V sin | in calculating the damping. This 
would occur for an azimuth angle of 270°, making the 
expression Qr — V 

Another possible approach is to use the value of | Qr + 
V sin y| corresponding to the average value over a 
cycle in the part of the blade azimuth under investiga- 
tion. Thus, for instance, if the second harmonic mo- 
ments were being investigated on the advancing side 
(¥ = Oto 180°) and a complete cycle of second harmonic 
moment spanned this range, |Qr + (2/r)V| would be 
used, 2/mr being the average value of sin y over this 
Similarly, for Y = 180° to 360°, |Qr — (2/r)V 


range. 
In general, if a cycle of second har- 


would be used. 
monic moment spanned the range of azimuth angle 


from yj, to } + 7, 


V viiter 
Qr + sin y dy 
Si 


would be used. For any harmonic, k, the formula be- 


comes 
kV vi + (24/k) 
Qr + 9 sin y dy 
us 


It should be noted that for higher harmonics, the num- 
ber of radians of azimuth per cycle, 27/k, becomes 
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TABLE 1 
SOLUTION FOR VALUES OF ©,' AND 0,1! FOR THE SECOND HARMONIC COSINE LOADIKG 
Section A 1 2 3 4 | 5 | 6 | > 
, : x2 x2@,! A. ete? Le x2e.! 
e s ead SS P] pe oe j 
El, eT, La oe ee ee uy" 
E 1 
“y J ere "9 
[(1)+(6)+(10)+(14) 2" 
+(15)],,, (xt0°") (VW) yx (2), (3), +(4),,, (4),x (5), (7),,,+¢4) - 
10 0 () ty) ty) 689.1 0) ty) 
9 102 4170 .425340 .425340 | 2079.6 884537064 0 
8 1316.537064 1200 1.579844 2.005184 | 3466.8 6951.571891 .425340 
7 9109.124779 476.9 4.344142 6.349326 | 4812.9 30558.671105 2.430524 
6 42540.370769 258.6 11.000940 17.350266 | 6159.3 106865.4934 8.779850 
5 160108.0476 163.5 26.177666 43.527932 | 7564.8 329280.099994 26.130116 
4 528767.8358 112.5 59.486382 | 103.014314 | 8722.0 898490. 8467 69.658048 
3 1545125.2166 82.38 |127.287415 | 230.301729 | 9560.2 2201730.5896 172 .672362 
2 4052658.4187 62.86 |254.750108 | 485.051837 |10139.6 4918231 .6064 402.974091 
1 9739687.1509 49.60 |483.088483 | 968.140320 |11092.9 | 10739483.7557 888.025928 
O | 24633513.5010 40.14 | 988.789232 |1956.929552 1856.166248 
8 9 10 "1 12 13 14 15 
; ; 
p , ie , 
202 2092 I 2092 I 20 Il =) II 20 II 
i 4x?70?m, | 4x7? m,U, 4x70 m,U; -2X?a, | U, -2X?Na,U, f  -2X?Na,U, oXS,.,° 
I*9 j=9 
7 * 
(7), x (8), (10),,,+(9), hak ve: | Gis, (14), ,, +13), 
10 0) ) 0 ft) 0 0 102 
9 618.0 ) 0 -584 ) ) () 330 
8 693.6 215.015824 295.015824 -596 ) 0 ) 546 
z 769.2 1869.559061 2164.574885 -590 0) ) ) 708 
6 897.6 7880.793360 10045 .368245 -564 .120895 -68.184780 -68.184780 725 
5 1124.4 29380.702430 39426 .070675 -518 . 894204 -4€3.197672 -531.3 485 
4 1157.2 80608 .293146 120034 .3638 -436 3.962035 1727.447269 58. 91 
3 1117.6 192978.6318 313012.9956 -328 13.614492 | -4465.553376 -6724.3831 | -486 
2 1158.8 466966 .3767 779979 .3723 -88 32.657539 | -3469.€63432 | -10214.2465 -968 
1 3813.2 | 3386220.4686 4166199.8409 i) 101.910242 ) -10214,.2465 -1643 
) 4166199.8409 233.735272 



































smaller, so that the approximation involved in averag- 
ing becomes more accurate. 

A further refinement is possible with some additional 
labor: This is the determination of an equivalent con- 
stant damping coefficient that would dissipate approxi- 
mately the same amount of energy per cycle as the 
actual damping. Here again, it would be necessary to 
select the part of the azimuth under consideration. 
A difficulty that arises in both this method of calcula- 
tion and the previous one is that, because of the phase 
shift introduced by aerodynamic damping (as will be 
shown later), the exact azimuth locations for the be- 
ginning and end of a cycle of kth harmonic blade vibra- 
tion are not known at the start of the computations. 
The location of the forcing moment cycle is known, 
however, and this may be used as a satisfactory ap- 
proximation in many cases. As a check on this, the 
damping factor may be recomputed on the basis of the 
calculated phase shifts, after the calculations based on 
the more approximate value of damping factor have 
been completed. It will usually be found that the 
change is small, so that nothing further will need to be 
done. In any case, too much concern over the accu- 
racy of this factor is not compatible with the order of ac- 
curacy of the other assumptions made in this analysis. 
All of the detailed calculations in this report are based 


on the assumption that a may be represented by its 
mean value over all azimuth positions—namely, 
a = Ao(p/2)cQr (4) 


Eq. (1) now becomes 


= (27 -) — Fir) ws +m = 
or? . or? of? 
p = +a = =f (5) 
or ot 


If to’ is zero, this equation is nothing more than the 
equation for the free vibration of the blade including 
A solution of Eq. (5) in closed form is not 
For this reason, it is desirable to 


air damping. 
at present possible. 
approach the problem from the energy standpoint com- 
monly used in vibration problems. 

Let the deflection of the blade be represented by the 
series 


x 


>» InPrn(1) (6) 
1 


Z_= 
n 
where g, is a function of time only and ¢, is a function 
of r only, which satisfies the boundary conditions for 
blade bending. The homogeneous differential equation 
for undamped blade vibration obtained from Eq. (5) by 
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ios TABLE 1 (CONTO) 
SOLUTION FOR VALUES OF e.! ano o,!! FOR THE SECOND HARMONIC COSINE LOADING 
section B 1 2 3 4 5 6 7 
El, EI, EI . i-% ET, : 
§*6 79 
[(1)+(6)+(10)+(14)),,,] (xt 0-8) (1), «(2),) (3), 404),,, (4), +(5), Cc ne 
10 0 b) vy) t) 689.1 t) vy) 
9 ) 4170 0 ) 2079.6 0 ty) 
8 ty) 1200 te) t) 3466.8 ) 0 
7 253502640 476.9 .120895 .120895 4812.9 581 .855546 0 
6 2522 .869986 258.6 .652414 773309 6159.3 4763 .042124 120895 
5 14033774662 163.5 2.294522 3.067831 7564.8 23207.527949 894204 
4 58530.0082 112.5 6.584626 9.652457 8722.0 84188.7300 3.962035 
3 198963 .2196 82.38 | 16.390590 | 26.043047 9560.2 248976.7379 13.614492 
2 576036.5299 02.86 | 36.209656 | 62.252703 | 10139.6 631217.5073 39.657539 
1 1416767.4856 49.60 | 70.271667 | 132.524370 | 11092.9 1470079 .5840 101 910242 
0 34849646532 40.14 | 139.886481 | 272.41085!1 234434612 
dam 
8 9 10 if 12 13 14 15 
: t 
1 | 4x?atm, | 4x?%m,u,!! axeat Y “ayu,"! +2x? Na, u,! 2x*N aU, ! ) 2xtaay,! 
78 j=9 
(906 
(7), (8), (10), 4,409), Re (11) x (12), (14) 54,4013), 
10 ft) v) 0 fv) 0) ty) ft) 
9 618.0 0 0 584 ty) 0) v) 
8 693.6 ) ) 596 425340 253 .502640 253.502640 
7 769.2 0 ) 590 2.430524 1434 ,009160 1687.511 800 
6 897.6 108.515352 108.515352 564 8.779850 4951. 835400 6639.347200 
5 | 1124.4 1005 .442978 1113.958330 518 26.130116 13535.400088 20174.747288 
4 | 1157.2 4584. 866902 5698. 825232 436 69.658048 30370.908928 50545.656216 
3 | 1117.6 | 15215.556259 20914,381491 328 172 .672362 56636 .534736 107182.1909 
2 1158.8 | 45955.156193 66869.5377 88 402.974091 35461 .720008 142643.9109 
1 3813.2 | 388604.1348 455473 .6725 ty) 888.025928 0 142643.9109 
fv) 455473 .6725 1856 .166248 
































setting fo’ and @ equal to zero has an infinitely de- 
numerable set of solutions,!! Z, = W,(r). A natural 
frequency, w,, corresponds to each solution Z,. The 
functions W,,(7) satisfy the orthogonality conditions, 


°R 
| mvV,V, dr = 0 
0 


The functions, V,,, are the usual vibration modes. 
modes may be used as the functions ¢,(7) in Eq. (6). 
Next, the loading, fo’, may be represented by a Fourier 


(when s * n) (7) 


These 


series of the type, 
to’ = to(r) + > [t.,(r) cos RQt + t,,(r) sin RQt] (8) 
k 1 


where ¢,,(r) and ¢,,(r) represent the kth harmonic 


cosine and sine blade loadings. 


When the loading is given by a series of the type rep- 
resented by Eq. (8), the principle of superposition may 
be used in solving the linear equation, Eq. (5). That 
is, the equation may be solved independently for each 
term in the loading series, and the solutions may be 
added to give the final result. The loading term, 4, 
merely corresponds to a static blade load for which the 
methods of solution are well known! so that it is not 


necessary to discuss this further. The kth harmonic 


cosine or sine blade loading may be represented by the 
real or imaginary part of a complex quantity of the type 
(10) 


be’ on Ww 


The solution of Eq. (5) with the above loading may be 


written as 


Z.= DY Ane — ™$,(r) (11) 
n=! 
and, accordingly, 
OZ,./dt = + ikOZ,\ 12) 
0°Z,./Ol? —k?0°Z,.4 “ 


Eq. (5) then becomes 


0°” 
or? 


EI — F(r) = te = + 
4 = r ) 4 
( or? \ or? f or 

Z(—k2Q2m + ikQa) = tie" (13) 


and, substituting Eq. (6) for Z, 


a °” , os) . 0*o,, 
’ EI — Fi(r) 
» 4 | or* ; or° 
p = + o,(—k2Q2m + ikQa) | = te (14) 
r 
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The principle of virtual work may now be applied. A displacement of the form ¢,(r) is assumed, and the virtual 
work is obtained by multiplying Eq. (14) by ¢, and integrating over the blade span, so that 


: TO? Ob, 0°>,, Wn ° ‘ 
YS dn f | (zr wi ) o, — Fi(r) * +p . o; tonbs(—R?Q?2m + ik) dr = f tip, ec dr (15) 
0 LOr? or? Or? or 0 


n=1 


It is now convenient to assume that the damping coupling terms 


R 
anita) f a(r)o,o; dr (n # s) 
0 





are negligibly small. If this is not true, then a system of simultaneous equations for the q,’s will result. The 
assumption that the damping coupling terms are negligible leads to a major simplification in the application of this 
analysis. In particular, it results in the possibility of expressing the amplification factor by a single relatively 
simple formula for the case where any particular forcing harmonic is near resonance with only one blade natural 
frequency. This is usually the case in typical helicopter blades because of the separation between blade natural 
frequencies. In this case, it may be shown, by means of a perturbation analysis similar to that given by Rayleigh," 
for a general damped dynamic system that the q’s other than the one corresponding to the mode near resonance 
will be considerably smaller in magnitude than the g corresponding to the mode near resonance, thus justifying 
the assumption that the damping coupling terms may be neglected. 
From the orthogonality condition, Eq. (7), 


R 
—k?Q? f mono; = 0 (2 ~ s) (16) 
0 


Also, since ¢,(7) is a solution of the homogeneous differential equation for blade vibration corresponding to the 


natural frequency w,, 


RTO? 26) 0*bn Od», . it e 
J [2 (zr rm — Fir) Sy? +p =] ¢; dr = —w,? . mo,o; dr = 0 (n x s) (17) 


Accordingly, all terms in the summation on the left-hand side of Eq. (15) vanish, except for those involving ¢, alone, 
) | g 


leading to the equation, 


RT O2/ 0%, 0°¢, Od; . ” 
qs [ | (er *) ¢; — Fir) _ o; + p . bs + o3°(—R?Q2? m + ita) | dr = f tid, e dr (18) 
JO Or or? or? or 0 


Integrating the first term in Eq. (18) by parts twice, 





R 9 9 9 R 9 R R 9 9 
Y /. a) ”) ( oe og, O°, f (=*) 
EI , dr = gy: El ~ -EI EI l (19 
I | or? _— or Or? J or Or? 6 ba 0 are)“ ) 
From the boundary conditions of the problem, 
o; = 0°,/dr? = 0 (at r = 0) ) 
0/dr [EI (0%,/ar?)]=Ol | (20) | 
0°4,/dr? = 0 peewee 


Thus, the first two terms on the right-hand side of Eq. (19) vanish so that 


R 9 9 R 9 9 
0° 0°, 0°, \* 
f (zr “6 dr = f El ( *) dr (21) 
0 Or or? 0 or? 
The second and third terms of Eq. (18) may be combined and integrated by parts to give 

R R “R 5 

O | 1) Obs O¢; | ; (o*): 
_ F dr = —¢,| F F 1 22 
£ ° | Fo | oar 6 | Fe) el + ; as (22) 


By consideration of the boundary condition, Eq. (20), and the additional condition F(R) = 0, it is found that the 


first term on the right-hand side of Eq. (22) vanishes and 


RO d¢, d¢,\° 
_ f | Fe S| o, dr = F(r) ( *.) dr (23) 
0 Or or 0 Or 
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Substituting Eqs. (21) and (23) into Eq. (18) 
*R 9 R -R 
qs | / EI ¢ a) dr + for o(% ") dr — k?Q? f mo,” dr + 1kQ I ads” ar| 
J 9 or? 0 0 


The natural frequency, w,, of the undamped rotating blade which corresponds to the normal mode, ¢,, is given by 


op >: 
| BI (> “)' rt fom (FS) dr 
0 or? 


Ws” — 
ia md,” dr 
0 


so that dividing Eq. ) through by 


f er (28) r+ fo ri) (2 ar 
0 or 0 


R 
io ad,” f tp er dr 
0 
gs} 1 — si R 2 R ° (26) 
0? 8 2 re) s zs 
var" mo,° dr f eI ( *) dr + f F(r) ) (*) dr 
0 or? 0 or 


R 
eikiat f tad, dr 
0 
“R 9 9 *R . ° 
0°; \* d¢,\? 
| EI ( *) dr +f Fir ( *) a 
ikQt — ids) ws Or? 0 or (an 
Lad 


ds — Ase = f 
ro ag,” dr 
7 tRQ 0 


1-4 — 


Ww,” Ws” r , 
md,” dr 
0 
-_ : 
;= J ag,” dr if md,” dr (28) 
0 / 0 
R *R ron \2 *R 2 
., (2¢s\" , (2s \" 
f ths ar/| f, pr (>) dr + / k (*) ar| 7 
(2% 


V [1 — (R?2?/w,?) ]? + (R??0,7/ws') 


ll 
= 
> 
= 
S 
. 
n 


results in 


Thus, 


and designating, 


there results 


SS 3 


6, = are tan [(kQa,/w,”)/1 — (R?2?/w,*)] (30) 


Returning to Eq. (11), the complete expression for the deflection of the blade under load, t,, is seen to be 


; i 
| t,o, dr 
0 
*R 0°6,, 2 *R : "6, " 
Z. = > / EI (S*) dr + / F(r) ( Ps ) ar ght ~- © 6/7) (31) 








If all the natural frequencies, w,, are extremely large compared to RQ, it can readily be seen that the denominators of 


all of the terms in the series approach unity and all the values of 6, approach zero. Eq. (31) then becomes 








820 JOURNAL OF THE AERONAUTICAL SCIENCES—DECEMBER, 1951 

























































































TABLE 2 | 
SOLUTION FOR VALUES OF f,1 amo f,'! FoR THE SECOND HARMONIC | 
Section A 1 2 3 4 5 6 7 
1 1 
x2 xz¢ 1 J ag 1 ea | 
. t," : f ott xt; Feet ve 9) fy v,? 
er ct, EI, ' EI, | 
ee j*9 aa 
((1)+(6)+(10)+(14)],,, (1) x (2), (4),,, +43), (4), (5), (7) 54,414), ,, | 
10 to) 0 1.0 689.1 -689.1 ty) 
9 -759.8 4170 -3.168366 4.168366 2079.6 -8668 .533934 -1.0 
8 -10117.033934 1200 -12.140441 -16.308807 3466 .8 -56539.371137 -5.168366 
7 -70755.911329 476.9 -33 743494 -50,052301 4812.9 -240896 .71 8679 +21.477173 
6 330985 .562742 258.6 -85 592867 -135.645167 6159.3 -835479 .279320 +71 .529.474 
5 -1244528 575977 163.5 -203 .480422 -339.125590 7564.8 256541 7.259450 -207.174 641 
4 -4103147.605490 112.5 -461 .604106 -800. 729695 8722.0 -6983964 .401098 -546 .300 230 
3 -11966382 .220528 82.38 -985.790567 | -1786.520262 9560.2 -17079491 .008772 -1347.029 926 
2 31334148 562866 62.86 | -1969.664578 | -3756.184840 10139.6 -38086211 . 803664 -3133.550188 
1 -75273759 .406854 49.60 | -3733.578466 | -7489.763306 11092.9 - 830831 95.377127 -6889.735 028 
0 -190482291 .433 40.14 -14379 .498 334 
8 9 10 "1 12 13 14 
1 ; 1 
2 I 202 I 2 II 2 II 2 II 
i | 4x?0%m, 4x70? mV, 4x? ) mV; -2X?Q a, vy -2X*NasVy } -2x?7Na,V; 
379 3*0 
, (Wy 
(7) x (8), (10), ,, +(9)4 aie oa (11) ,x (12), (14), ,,#(13), 
10 70.7 0 0 0 
9 618.0 -618.0 -688.7 -584 0 0 0 
8 693.6 -3584.778658 -4273 .478658 +596 -.291900 173 .972400 173 .972400 
7 769.2 -16520.241256 -20793.719914 -590 2.181042 1286.814780 1460. 787180 
6 897.6 -64204 855461 -84998 575375 -564 -9.705770 5474 .054280 6934 .841460 
5 1124.4 -232947.166228 -317945.741603 -518 «34 380560 17809. 130080 24743 .971540 
4 1157.2 -632178.626677 -950124 368280 -436 -105.757300 46110.182800 70854 . 154340 
3 1117.6 -1505440.645298 +2455565 .013578 -328 -294.010749 96435 .525672 167289.680012 
2 1158.8 | --3631157.957854 -6086722 .971432 -88 -750.389217 66034 .251096 233323 .931108 
1 3813.2 | -26271937.608770 -32358660 .580202 0 -1770.125241 0 233323 ..931108 
R 
: td, dr 
0 
Pere, ikQt 90) 
Z,' = > Rk /d%,\2 R oe? |e. bal”) (32) 
. n </ n 
mg” EI > dr + F(r) ( dr 
0 2 0 or 


This is the result that would be obtained from an energy solution! for blade bending under a periodic loading, 


( . . rr . . . . oe . 
t,e"*, if dynamic effects were neglected. Thus, it is seen that consideration of the dynamic effects leads to an | 
amplification of each term of the series by the factor 
F, = 1/V [1 — (R?Q?/wn?)]? + (R72? on?/an*) (33) 
and a lag in phase between the bending moment and the periodic aerodynamic loading of amount | 
6, = arc tan (RkQa,/w,”)/[1 — (R?2?/w,”) | (34) | 


Since the set of orthogonal functions, ¢,(7), corre- nitude, since they are excited by higher harmonics of 
sponding to a number of the lowest blade vibration the loading which theoretically are extremely small 
modes cannot be easily obtained in most cases, it is (see, however, reference 15). Therefore, the “‘static’’ 
usually not practical to solve for the deflection, Z,’, in moments corresponding to Eq. (32) can be determined 
the form of the series Eq. (32). However, it is often by a tabular method similar to that used in reference 
found that, for a typical rotor blade, the first bending 3, and the dynamic effects can then be considered as 
mode (n = 1) is predominant in blade bending and the amplifying these ‘‘static’’ moments by the factor of 
higher modes contribute moments of negligible mag- Eq. (33) with m = 1. The corresponding phase lag 
















































































DYNAMIC EFFECTS OF ROTOR BLADE BENDING 821 
—_—_ TABLE 2 (cCoONTD) 
SOLUTION FOR VALUES OF #,' AND f#,'' FOR THE SECOND HARMONIC 
Section B 1 2 3 4 5 6 7 
1 1 
x2 x? I! x2¢ tt x2¢ tt 
i €** gad —+_ —- eee Fi. rina vt 
EI El, 729 El, a" El, 
((1)+(6)+(10)+(14)],,, (ad ete), 1O),,, 219), (4), x (5), (Hien 
10 0 0 te) 0 689.1 0 te) 
a -70.0 4170 -.291900 -.291900| 2079.6 -607.035240 0 
8 -1331.035240 1200 -1.597242 -1.889142 | 3466.8 -6549.278484 -.291900 
7 -11817.121700 476.9 -5.635585 -7.524728| 4812.9 -36215.761596 -2.181042 
6 -66318.8806 83 258.6 -17.150062 -24.674790| 6159.3 -151979.435106 -9.705770 
5 -285638.835077 163.5 -46.701950 -71.376740| 7564.8 -539950.760543 -34.380560 
4 -1038904.080558 112.5 -116.876709| -188.253449| 8722.0 -1641946.580259 -105.757300 
3 +3254734.393791 82.38) -268.125019| -456.378468| 9560.2 -4363069.431112 -294.010749 
2 -8962099.203231 62.86| -563.357556/-1019.736024 |10339.6 | -10339715.388950 -750.389217 
1 -21791413.995169 49.60 |-1080.854134/-2100.590158/11092.9 | -23301636.563678 -1770.125241 
0 -54332491 .530800 40.14 -3870.715399 
8 9 10 11 12 13 14 
1 1 
i 4x70? m, 4x?n?m,v,t! 4x22 ) a,v,"! 2x7 a, ¥,* 2x?0a,V,! ) 2x20 a,V;! 
j=*9 3*9 
1 (70, = 
(7) x (8), (10),,,4(9), Pam a (11) ,x (12), (14),,, #013), 
10 fe) 0 fe) 0 0 -70.0 
9 618.0 fe) 0 584 -1.0 -584.0 -654.0 
8 693.6 +202 .461840 +202 .461840 596 -5.160366 -3080.346136 +3734.346136 
7 769.2 -1677.657506 -1880.119346 590 -21.477173 12671 ..531905 -16405.878041 
6 897.6 -8711.898817 -10592 ,018163 564 -71.529474 -40342 .623084 -56748.501125 
5 1124.4 -38657.501664 -49249.519827 518 -207.174641 | -107316.463986 -164064.965111 
4 1157.2 -122382.347560 -171631 .867387 436 -546.300230| -238186.900476 - 402251 .865587 
3 1117.6 -328586 .413082 -500218.280469 328 -1347.029926| -441825.815728 -844077.681315 
2 1158.8 -869551.024660 | -1369769.305129 88 -3133.550188| -275752.416544 | -1119830.097859 
1 3813.2 -6749841 .568981 | -8119610.874110 0 -6889.735028 0 -1119830.097859 





























can likewise be approximated by use of Eq. (34), set- 


ting 2 = 1. 
(o, = QO), 


simpler factor used by Johnson and Mayne. 


It should be noted that, for no damping 
the amplification factor, F;., reduces to the 


This 


same procedure can, of course, be applied to resonance 


with any of the first few bending modes. 


A more accurate method of applying the amplification 


factor may be employed based on Eq. 


(31). This 


method is particularly appropriate when the vibration 


mode under consideration, ¢), is known from a detailed 


vibration analysis such as Myklestad’s tabular analy- 


sis'* or a Rayleigh-Ritz solution. 


If the forcing har- 


monics being considered are near resonance only with 
w, the amplification factor and phase shift should, 
strictly speaking, be applied only to the ¢; component 
of the static deflection, Z,’, which, according to Eq. 


(32), is given by the series, 


tRQU 


Z.’ = DY and,(rye 
1 


n” 


where 


“R 
| Ln dr 
0 


(35) 


(36) 


an, = R a ha R : 
ra) n ra) n % 
[ rt ( *) dr + F(r) ( $s) dr 
/0 or? 0 or 


Using Eq. (25), 


mR PR 
a, = J typ, dr/w," J md, dr (37) 
0 0 
An equivalent formula, 
R °R 
any = mZ, cy dr | md,° dr (38) 
0 


in terms of Z,’ instead of ¢,, which is sometimes more 
convenient, may be obtained by multiplying Eq. (35) 
through by mq, integrating over the blade, and using 
the orthogonality conditions. Then, *Z,’, the static 
deflection without the first mode component, is given 


by 


“Z7’ = Zz. > ay,o,e (39) 


The deflection in the dynamic case is then given by 


i (RQ 5) 


Z, = *Z,' + ay,F die (40) 


where F; and 6 are given by Eqs. (33) and (34), respec- 
tively, using m = 1. The static blade bending moment 
less the first mode component is given by 


One ee 
*M' = El Y= M’ — a, EI 


(41) 
or? or? 
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Fic. 1. Natural blade bending frequency and harmonic load- 


ing frequencies vs. rotor angular velocity for a typical 19-ft. rotor 
blade. 


The blade bending moment in the dynamic case is given 


by 


0? — 
M = *M' + a,,F;, (El) ™  -* (42) 
or? 
This may be written 
M = *M'+ M'Fe-™” (43) 


where ./,’ is the first mode component of the static 


bending moment—namely, 


0°) ikQt 
é 
or? 


M,’ = Aa, (El) (44) 


When ¢; has been calculated by approximate or tabular 
methods, as is usually the case, double differentiation 
to obtain 0°¢,/Or? will greatly diminish the accuracy 
of the results. This may be avoided by taking ad- 
vantage of the fact that ¢;, being a vibration mode of 
the rotating blade, satisfies the equation 


0° 0’ #) 0°¢1 Od: 
E — a i 
= ( yf Dp? F(r) Sy? + p Sp ¢ = 0 


(45) 


Thus, 
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7o o 0° 6 Op 
" SS * ri L rns ne r net 
R (R 
f i mo"; dr; drz (46) 


where 7; and 72 are merely running variables and dis- 
appear after integration. Making use of the fact that 


— |OF(r)/Or] and integrating by parts, this becomes 


Oo” R 
EI =— 7 P(r) [di(r1) — o1(7)] dry + 


or? 
R 
f m oa" oi(%1, — 7) dry (47) 


This equation provides the required value of E/(0°¢,+ 
Or’) for substitution into Eq. (41) in terms of integra- 
tions of ¢@ only. The practical application of the 
method under discussion here is carried out by using 
Eqs. (41), (48), (44), (47), and either (37) or (38). 


Application 


For typical rotor blades and rotor r.p.m.’s, it is ob- 
served, as previously stated, that the fundamental 
bending frequency, w, is of the same order of magni- 
tude as the frequencies of the lower harmonics of the 
loading. Therefore, for an example of the preceding 
theory, it is assumed that the first bending mode is 
predominant in the blade deformations and the higher 
modes contribute negligible effects. Fig. 1 shows the 
first natural bending frequency and also the important 
aerodynamic loading frequencies for various rotor speeds 
as calculated for a typical blade. A diagram of this 
type indicates the forcing harmonics that are near 
resonance and that should be made for any case under 
consideration. 

Based on the above conditions, the detailed proce- 
dure for applying the theoretical results to a practical 
problem is as follows: 

Step Calculate the 
moments by a tabular method. 

calculation is indicated in reference 3, page 2 

Step 2.—Calculate the blade’s first aa bending 
frequency, w. Reference 3, page 20, evaluates w, for 
a typical blade. 

Step 3—Calculate the value of o; by use of Eq. (25) 
and an assumed deflection mode. 

Step 4.—Calculate the amplification factors for the 
various harmonics by use of Eq. (33) with m = 1. 

Step 5.—Calculate the phase shifts for the various 
harmonics by use of Eq. (34) with” = 1. 

Step 6.—Obtain the dynamic moments by applying 
the amplification factors and phase shifts to the tabu- 


various harmonic ‘‘static’’ 


The — for this 


lar-‘‘static’’ moments. tT 

+ If the mode, ¢;, has been calculated rather than assumed, it 
will often be worth while to carry out this step by application of 
Eqs. (41), (43), (44), (47), and either (37) or (38). This will be 
particularly true if the dynamic effects are found to be critical 


in the stress analysis. 
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To evaluate the effects of aerodynamic damping on 
the bending moments, the amplification factors and 
phase lags were ascertained for the typical rotor blade 
in a forward flight condition (120 m.p.h., 252 r.p.m.). 
This is the same blade and flight condition treated by 
Johnson and Mayne* in which blade damping was 
neglected. For the first harmonic, an amplification 
factor of 1.184 was calculated as compared to 1.217 
obtained in reference 3. For the second harmonic, 
which as shown in Fig. 1 is the closest to resonance, a 
factor of 2.02 was obtained as compared to 3.472 pre- 
sented by reference 3. Moreover, consideration of the 
damping effect yields a phase lag of 14° for the first 
harmonic and 54° for the second harmonic. Reference 
3 indicates zero phase lag, since as stated previously, no 
damping is considered. Thus, as would be expected, 
the significant effects of aerodynamic damping are to 
decrease the amplification of the bending moments and 
also to introduce a phase angle between the bending 
moment and the periodic aerodynamic loadings. 

It is interesting to note here that, as shown in Fig. 1, 
resonance would occur for a rotor angular speed of Q = 
11.7 rad. per sec. The maximum amplification of the 
bending moments would appear to occur at this point. 
According to the simple undamped single degree of 
freedom analysis,* the amplification factor for this 
speed would be infinite. However, if the aerodynamic 
damping is included in the manner suggested, this in- 
finite factor is reduced to a value of approximately 


TABULAR METHOD INCLUDING AERODYNAMIC 


DAMPING 


THE 


The more exact tabular method devised by Johnson 
and Mayne’ has been extended to include the important 
aerodynamic damping terms. The solution that fol- 
lows is applicable to any hinged rotor blade. 


Theoretical Development 


The net normal force distribution acting on a rigid 
blade hinged to the rotor hub is a result of the inertia 
and aerodynamic forces on the blade. The inertia 
forces are functions of the blade’s mass distribution, 
angular velocity, and flapping motions. The various 
harmonic flapping motions are readily determined for 
any flight condition.'‘ Also, in conjunction with the 
evaluation of the flapping motions, the various harmonic 
aerodynamic loads can be ascertained. Consequently, 
the net normal force distribution for the various har- 
monics that act transversely on the blade may be de- 
termined. These net normal force distributions are 
considered as the forcing functions that tend to deflect 
the blade. Naturally, the elastic blade on deflecting 
introduces additional inertia and aerodynamic forces. 

The basic equation governing the elastic deflection of 


a rotor blade is 


ROTOR 


BLADE BENDING 823 


dM/dr = H(r) (48) 


The total shear force, 7, at any station, r, of an elastic 
blade is considered to consist of the following: 

(1) S(r) = the net shear acting on the rigid blade 
specified above. 

(2) F(r) (dz/dr) = shear component of the centrif- 
ugal force introduced by bending of the blade. 

(3) JZ (r) = shear of the inertia forces introduced 
by blade bending. 

(4) D(r) = shear of the aerodynamic forces due to 
blade bending. This force has the characteristic of a 
damping force, being proportional to the deflection 
velocity. 

Accordingly, 

H = $+ F(dZ/dr) +T+D (49) 
Substituting Eq. (49) into (48), there results 

dM/dr = S + F(dZ/dr) + T+ D (50) 

The kth harmonic deflection of the blade is of the form 


Z = Z' cos kQt + Z" sin ROL (51) 


rz1 ° — 
where Z' and Z" are functions only of r. Using the 


above relations, the various shears take the form 
dz dz! dz" 
F( )- F cos RQt + - sin RQt 
dr dr dr 


R ss 
T = f mZ dr = 


R 
—wne f m([Z' cos kQt + Z" sin RQ] dr 
r 


R * 
p- f aL dr = 


R 
io f a| —Z' sin kot +Z" cos RQt] dr 
S = S' cos kt + S" sin kM ) 
The kth harmonic bending moments can therefore be 
expressed in the form 
M = M' cos kOt + M" sin k0e (53) 


where M' and M" are functions only of r. Substitut- 
ing Eqs. (52) and (53) into Eq. (50) and separating the 


sine and cosine terms gives 


dM' _sl4iF dZ' a 
dr dr 
R *R 
eo f mZ' dr +kQ / aZ" dr 
a (54) 

I zil 
= = S'4F- = — RQ? X 
dr dr 


R "R 
" onl 
z mZ" dr — kQ | aZ dr ; 
’ r 
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TABLE 3 
SOLUTION FOR VALUES OF g,! AMD g, 1! FOR THE SECOND HARMONIC 
1 2 3 4 5 6 7 
Section A 
1 : 2 I - 2 1 
I x? x2 x*g x2g 
’ at Es, meee se "9 Fy? — wt 
et, ET Wee c- 8 
COV#(G)H(VO)H(14)],,, | CxFO-F DT CTD (2d, | (4d, F038), (4) x (5), Fg st Deas 
10 0) ) 0) ty) 689.1 ft) 1.0 
9 - 282.8 4170 1.179276 1.179276 2079.6 2452 .422370 1.0 
8 3636 .022369 1200 4.363227 5.542503 3466.8 19214.748856 2.179276 
7 24567.227459 476.9 11.716111 17.258614 4812.9 83063. 981482 7.721779 
6 111380.413751 258.6 28.802975 46 .061589 6159.3] 283707.142750 24.980392 
5 407602 .104060 163.5 66.642944 | 112.704533 7564.8] 852587.248212 71.04198 
4 1314484.515292 112.5 | 147.879508 | 260.584041 8722.0] 2272814.002113 183.746514 
3 3766781 .395919 82.38] 310.307451 | 570.891492 9560.2] 5457836.841723 444 .330554 
2 9734062 .109354 62.86] 611.883144 | 1182.774636 | 10139.6 |11992861 .699186 1015.222046 
1 23306502 .156872 49.60|1156.002507 | 2338.777143 | 11092.9 |25943820.969585 2197.996682 
O | 59211302.422591 40.14 4536 .773825 
8 9 10 "1 12 13 14 
H i 
i | 4x?0?m, 4x70? mw, ! axa) m,W,! -2x?N a, w,** -2x*M awit p -2x?2a,witt 
3=9 j=9 
(7), 
(7), (8), (10), 4, +09), aitiien asi (Pt) nf12%,, (14), ,, +113), 
10 282.8 ty) ty) 0 
9 618.0 618.0 900.8 -584 0 0 0 
8 693.6 1511.545834 2412.345834 -596 1.167600 -695 .889600 -695.889600 
7 769.2 5939.592286 8351 .938120 -590 6.621769] -3906.843710 -4602.733310 
6 897.6 22422 .400273 30744 .338393 -564 24.214641 | -13657.057524 -18259.790834 
5 | 1124.4 79879.603521 110653.941914 -518 73.550170| -38098.988060 -56358.778894 
4 | 1157.2 | 212631.465654 323285.407568 -436 200.559C60| -87443.750160 - 143802 .529054 
3] 1117.6 | 496583.827486 819869.235054 -328 507 .996446 |-166622 .834288 -310425.363342 
2 | 1158.8 | 1176439.307252 | 1996308.542306 +88 | 1208.009439 |-106304 .830632 -416730.193974 
1 3813.2 | 8381400.947802 | 10377709.490108 0 2703 .283853 0 -416730.193974 
































equations. Dividing the blade into an integral number 
of parts, m, each of length X, the derivative dM//dr 
at the midpoint of an interval may be expressed approxi- 


Inspection of Eq. (54) shows that the damping effect 
due to the sine deflection enters into the determination 
of the cosine moments and, similarly, that the damping 
effect due to cosine deflection enters into the determi- mately by 
nation of sine moments. Thus, it is to be expected that CNA net, = ME OD 
the solutions of the above equations must be obtained P 


by means of a simultaneous treatment of the two Substituting Eq. (55) into Eq. (54) gives 


I I ro I - 7 dZ' 
M; — = M, = X S; ana = F; — (1/2) X + 
s — (1/2) 


dr 
R *R 
Xk202 mZ" dr — XkQ | aZ" ay 
r=1 (1/2 r =4— (1/2) 
sa (56) 
éa™ 
M; - r ‘ii M" ~AnN-~«@ 2) —F;-an2% dr ) ? 


(1/2) 
> 


4 R 
- 1902 r , I 
X k?Q? / mZ" dr + XkQ / aZ dr 
af? t (1/2) ’ # — (1/2) 


The integrations contained in Eqs. (56) may be approximated by summations. It is seen that the values of 7 


and D are desired at the half stations. These values are obtained by determining the values of k?Q?mZ and 
kQaZ at the unit stations, at which point they represent the mean value between the half stations, and by multi- 
However, the effect of the portion of the blade at the tip be- 


plying these mean values by the interval length X. 
For this portion, the values of k°Q?mZ 


tween station m — (1/2) and n must also be included in the summation. 

















er 
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TABLE 3 (conto) 
SOLUTION FOR VALUES OF o,* ano ~-- FOR THE SECOND HARMONSC 



























































Section 8B ' 2 3 4 5 6 7 
i 1 
at! xt =. a = a oe) ee wt 
1 Sig 4 Ely : ; ‘, 
a*e j*9 
[(1)+(6)+(10)+(14)],,,] (X10- 8) (30,0120, } (43,,, +193, (4),x (5), CF cig hi 
10 0 0 0 689.1 r) 0 
9 280.0 4170 1.167600 1.167600 | 2079.6 2428.140960 0 
e 3572.140960 1200 4.286569 5.454169 | 3466.8 18908.513089 1.167600 
7 25453 .349905 476.9 12138703 17.592872 4812.9 84672.733649 6.621769 
6 122748 .093642 258.6 31.742657 49.339529 | 6159.3 303872 .323770 24.214641 
5 475066 .430609 163.5 77.6733€1 | 127.008890 7564.8 960796 .851072 73 .550170 
4 1603808. 852223 112.5 | 180.428496] 307.437386 8722.0 | 2681468.880692 200.559060 
3 4765423 .727662 82.38] 392.575607| 700.012993 9560.2 | 6692264.215679 507. 996446 
2 | 12651311.187948 62.86 | 795.261421 | 1495.274414 | 10139.6 | 15161484, 448194 1208.009439 
1 | 30495599.758736 49.60 |1512.581748 | 3007.856162 | 11092.9 | 33365847.619450 | 2703.283853 
0 | 76852413.489040 40.14 5711.140015 
8 9 10 "1 12 13 14 
1 i 
1 | 4x2%m, | 4x?Q?%m,w,?! arer om," 2x?N ay w,! 2x?Na,wt 7 ney 
j=9 ites 
(7),« (8), (10), ,, ¢09), rrem (11) x (12), (14),,, 4013), 
TABLE 4A 
10 1.0 280.0 
9 618.0 C) t) 504 1.0 584.0 864.0 
8 693.6 809.847360 809 .847360 596 2.179276 1298. 848496 2162.848496 
7 769.2 5093 464715 5903.312075 590 7.721779 4555.849517|  6718.698013 
6 897.6 21735.061 762 27638.373837 564 24 .980392 14088.941347] 20807.639360 
5 1124.4 82699.811148] 110338.184985 518 71.041981 36799.746197| 57607.385557 
4 1157.2 232086944232] 342425.129217 436 183.746514 80113.479973] 137720.865530 
3 1117.6 567736.828050| 910161.957267 328 444330554 145740.421810| 283461 .287340 
2 1158.8 | 1399841 .337913] 2310003.295180 88 |1015.222046 89339.540074| 372800.827414 
1 3613.2 | 10308161 .988260| 1261 8165.283440 oO 2197.996682 0 372800.827414 
0 


























and kQaZ are evaluated at the midpoint of the interval [at m — (1/4)]|, and the values so obtained are multiplied 


by the half interval length X/2. Accordingly, 7 and D may be written in the following form: 


i 


>R y 
T = Rk? / mZ, dr = k*Q? 5 Mn — 1/4) Zn -/sy FRAX YD mdz; 
Jr=i (1/2) 2 j n—1 | 
: (07) 
" ¥ i 
D = kQ / aZ dr = kQ > an — ays) Sn — cay + ROX > ad; 
Jt $ — (1/2) - jzn 1 


Substituting Eq. (57) into Eq. (56), the following expressions are obtained : 


y? 
9, ” zi y I 
) + k°0? > Mn -( Zn - ass) + | 
# — (1/2) = 


; ' _ : , dZ' 
M = M, = = S; (1/2 = F; (1/2) X 
dr 


1 €. i 


9”)9 Y2 zl a r ro mm | 
R2Q2X <j y m,Z; —_ RQ 9 Qn — (1/4) Ae = isa RQX ‘J > a;Z; 
J n—l ~ J n—1 rt (58) 


rs | 


ail 

dZ ae : - 
) is k°0? 9 My - (1 yZn — (1/4) + 
t — (1/2) cd 


dr 


M;. = M/" — XS; — ayn" — Fi ~ ayn X 


1 


i X2 

9-9 \Y79 7 ll = Y I r? 7 I 

R°O2X? DD mZy + RQ on — ayn - (1/4) + ROX? > ayy 
j=n-1 «= j=n-1 


In order that Eqs. (58) may be written only in terms of moment parameters and the boundary conditions, the fol- 
lowing relationships are used :* 
; 


(dZ/dr); — 2) = (€Z/dr)n — 72) — d [MjX/(ED,) 


j=n-l1 


(59) 
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and 


Z;-1 = Z, — X (dZ/dr)n — ayy + YX (MXV(ED,) (60) 
a—] 


j 


To simplify the notation, let 





y' = X (dZ"/dr), — 12); yt = X dZ"/dr)» — (12) (61) 
Thus, it follows that 
Zn — (1 4)" a z; — (y' 4); Zn qd . = FR aed (yl! 4) (62) 
Substituting Eqs. (59), (61), and (62) into Eqs. (58) gives 
M;_,) = M/ — XS > me + #0? m zi —%| 4 
==] et toed | (1/2) (El), 2 n — (1/4) “n 4 
j he 
5 . X? yi 5 
R?02X m,Z,;\ — kQ 5 Om - (1/4) Z.° 4 — RQX? a Yt 
j=n-1 ss j=n-1 (63) 
xu x? yl 
M.S ol +h. aa” - Faas = > —— | 4 £2? — mg | Z,1 — — [+ 
(El); 2 4 
j=n-1 = 
i | 1 
202 ¥2 me a= call Kg or? = 
R°02X? my + RQ om — 1/4) | Zn — 4 | + RQX? a2; 
j=n-1 . j=n-l 
Also, by substituting Eqs. (61) into Eq. (60), it follows that 
Z-{=Zi-—y+ DY (M)'X2/(ED);) 
j n—1 
(64) 
Z- =Z/'-y"+ DY (M;"X)/(ED, 
j=n—1 
It is noticed that Eqs. (63), together with Eqs. (64), express the moments at any station (7 — 1) in terms of the 


: ; HI og! - ti . = ae eee , 
moments at the station (7) and the terms y', 2 ge Sc ee As in reference 3, it is found expedient to repre- 


sent the moments and deflection by the four oneages —", 
I I I7 I a yz I 
M, ms o£ + fi'y roy £i Ze = y' = Be Se | 
I I AL 1 Iu I7 I 
M," =e)" + fily' + 902, + fey + gi Zp 


Z;} = U? 4. is + Ww,'z,! ae yilyt aa WUZ II (65) 
; Bas —_, > +. V ity! 7s fy 4. V iy" +4. WIZ," | 
Substituting Eqs. (65) into Eqs. (63) and equating the coefficients of y', yl, Z,', and Z,,"' gives 
I I 3 ~ ws’ —_ I stain + rl °F 7 II 
6-1 =4 + Fa YD = —XS;~ ayy + ROX? DY mU; —kOX? DY al (66a) 
J (EI), j=n-1 jz=un-l 
—_— ot : : X°f;! — My, — (1/4) 2 1 ae = 
fier =fi + Fe -ayl(-1+ > + R?02X? : + Do mV; ) — ROX? SS al 
jon—1 (EI); 5 j=n—1 j=n-1 
(66b) 
I I . X°*g;' Oe iy iain: ae 1 ro , , il 
Qi-1. =e +Fi-aa >» Stee > Mn — (1/4) + R°O2X2? SO mW, —kOX? YO al (66c) 
Janne GED e - j=n-!l j=n-1 
1 Ye II i 1 
e-visel+ Fa SY = +t khoex? Y mu" + ROX? DY aU — XS; — 2)" (67a) 
j=n-—1 (EI); j=n-—1 j=n-1 
; II me! , : im 202 V2 Il r9 an 1/4) >I i 
fier’ =fi' + F-a DY = + kX? DY mV," + OX? | - + >> aV;') (67b) 
j=n-1 (EI); j=n-1 5 j=n-—-l 


a ene i x 
g-1 =a +R ay Dd aie + k202X? SY mW," + ko > On — (1/4) + ROX? YS aW;! (67¢) 
n—1 “Lt )j j=n-1 “= j=n-l1 





(60) 


(61) 


(62) 


(63) 


64) 


the 
re- 


a) 


b) 
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TABLE 4 
CALCULATION OF MOMENTS DUE TO THE SECOND HARMONIC COSINE LOADING 

section A i e,! f,! y! o,! ‘7 -f,'! y"! -9,!! z,1! m,! 

10 0 0 0 0 0 8] 

9 102.0 -6.909458 109.429778 32 .039894 -251.518549 -20.958335 

8 1316.537064 -92.002131 1406 .962943 604.231835 -3285 .330562 -44.60085!1 

7 9109.124797 +643 .439045 9506.316279 5408 .847578 -23409.677635 -28.828044 

6 42540.370769 -3009.911546 43098.776292 30354 .998979 -112892 .539219 +91 .695275 

5 160108.0476 -11317.4753 157722.092311 130740.5441 +436922.9212 330.2875 

4 528767.8358 -37313.1421 508641 .2606 475519.6709 -1475037.6023 578.0229 

3 1545125.2166 -108819.7071 1457560.2950 1489733 .5146 -4382803 .5871 795.7320 

2 4052658.4187 -284946 .0102 3766606 .2742 4102067.3050 -11635526.0771 859.9106 

1 9739687.1509 -684523 .3842 901 8477 .2062 9974208 .5926 -28047080.7303 768.8352 

0 24633513.5010 |-1732205.0046 22911880.0261 24868675.5256 - 70681 864.3502 -.3021 
Section B i e,U ¢ ti yl git Zz f,? yl of 27 ae 

10 0 0 0 0 1?) .¢) 

9 i?) -.636565 108.346315 -347.770167 260.093735 20.033318 

8 0 -12.104148 1382 .243960 -4630.695687 3344 .082875 83 .527000 

7 253.502640 -107.462364 9849.202355 -32385.884593 22594.702754 204 .060792 

6 2522 .869986 -603 .089642 47497 .512803 -151496.320738 102437.580533 358.552942 

5 14033 .774662 -2597.538153 183827.4893 -569636.6293 374875.3659 502.4624 

4 58530.0082 -9447.5703 620595 .6380 -1878063 .0806 1208943 .3755 558.3708 

3 198963 .2196 -29597.8548 1843986 .0679 -5477166.0251 3464343 .1427 528.5503 

2 576036 .5299 -81499,.4033 4895439.0843 -14342040.1224 8952505 .5405 441.6290 

1 1416767.4858 -198166.4337 11800306.6038 -34453761 .3781 21435202 .2162 348.4940 

0] 3484964 .6532 -494087.9965 29738127.7817 -87186178.3893 54457173.8981 -.0528 














T = +,009093785 y'? = +,457712776 


y 
z= +.386951124 


z,21 = +.919709104 





Similarly, substituting Eqs. (65) into Eqs. (64) gives 


1 a 2 I 
U;- "s = U;' + > . 


: (68a) 
j=n-1 (EI); 


—— - . xy; :; 
VV; =Vi + ~t+ & (68b) 
J 


~ (an, 
W,_{ = Wi + > | a (68e) 
U;_ = ue + pos ae (69a) 
Vj; = V+ E | ae (69b) 
WwW," =we t+ oe (69e) 


It is now noticed that by a step-by-step procedure, 
which may be set up in a tabular form, the values of 
e, and e,'' can be obtained by use of Eqs. (66a), (67a), 
(68a), (69a), and the following boundary conditions: 


l Py - l the = en aad e, | aad 0 


(since ati = n; M,! = M,"! =0; Z,'=Z,' and Z,"' = 
Z,/!). The values of f;' and f;'' may be similarly ob- 
tained by use of Eqs. (66b), (67b), (68b), and (69b). 
The boundary conditions are, at 7 = n, 


ry," _ = a 4, sie tn = 0 


Similarly, the values of g;' and g;'' may be obtained by 
use of Eqs. (66c), (67c), (68c), and (69c). The cor- 
responding boundary conditions are, at 7 = n, 


W, -” i, Ww," _ Ln = a. = 0 


After obtaining the values of the various coefficients, 
the parameters } ee Aad y', and yl can be evaluated 
by use of the boundary conditions that the deflections 
and moments at the hinged end must be equal to zero. 


Accordingly, Eq. (65) yields 


0 = ey! + foly' + go'Z,' — forty — go'Z," 

0 ei! + folly! + gz, 1 + fly! + glZ," | 

0 Us + Very! + We'Z,! — Volty"™ — wUz," 

0 = ye" + Voity! 4. fg & + Voly™! +. W,'Z,"' 
(70) 


I 


The values Z,,', Z,,"', y', and yl are then used, together 
with the coefficients in Eq. (65), to calculate the moment 
at any station, 7, of the blade. 


Application 

In order to illustrate the application of the tabular 
method, as well as to check the accuracy of the more 
approximate method, the typical blade was analyzed 
by this method. The tabular procedure for deter- 
mining the moments due to the different harmonic 
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: “alne f, > glgl Jt iT 
od Es Step 4. Solve for the Parameters Z,,Z,, y', and y 
Tt (ates eae Having obtained the various coefficients (e, f, g, U, V, 
/ satin ‘ r\ P P 
j ; and |W), the parameters a, ==. y,! and yl may be 
“| Taaucan ra determined by simultaneous solution of Eq. (70). For 
30 AMPLIFICATION FACTOR. | rat . 
{ manne Ls the present example, the solution yields parameters 
oe J having the values of 


DEGREE 
= 
‘ 























“ 
< 
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+ 
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x z TABULAR DYNAMIC 
ga AMPLIFICATION FACTOR? 

= (MEGLECT AERODYNAMIC 
> DAM PING ) \ 
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N IN } 30) \ 

TABULA : | hs , A 

AADLitica TION FACTOR SS = 

cape en oaunn = 

° , if 2 3 r) 3 * 7 r) ¥ 
SPANWISE STATION (x } 


SPANWISE. STATION ({) 


Maximum station moments and phase angles due to 
rotor blade. For- 


Fic. 2. 
second harmonic loading for a typical 19-ft. 
ward flight, 120 m.p.h., 252 r.p.m. 


loadings is identical. Independent solutions were ob- 
tained for the effect of the second harmonic (k = 2) co- 
sine and sine components, and the former is presented 
herein. However, it should be noticed that both load 
components of a harmonic may be treated simultane- 
ously in the same tables. Consequently, as is shown, 
four tables and a check table are necessary for each har- 
monic under investigation. The detail procedure is as 
follows: 

Step Compute Table 1.—This table is set up to 
satisfy Eqs. (66a), (67a), (68a), and (69a). As pre- 
viously stated in the example second herein, S" is 
not considered. The blade properties, X?/(E/);, 
F;_ (12), R°?X?Q?my, RX°Qa;, and XS;_(,/2)' were evalu- 
ated and introduced in the appropriate columns. The 
starting conditions are ey! = ey"! Uw! = U,y'! = 0. 
The procedure then followed is indicated at the heading 
of each column, The two sections (A and B) of the 
table must be solved simultaneously as required in col- 
umn (12). 

Step 2. Compute Table 2.—This table is set up so as 
to satisfy Eqs. (66b), (67b), (68b), and (69b). It 
yields values of fy iP. Vi, and V,''. The starting 
conditions are Vy»! = V,!! = fio! = te" = & fe 
should be noticed that in Section A, the —70.7 
[= —(1/8) k?Q?X?m,,_(1/4)] is introduced in column 
(10), rowz = 10, in order to conform to Eq. (66b). 
Similarly, in Section B, the value —70 [= —(1/8) X 
ROX? a, — (4 row2 = 10, 
to conform to Eq. (67b). 

Step 3. Compute Table 3. 
satisfy Eqs. (66c), (67c), (68c), and (69)c. 


value 


,)] is introduced in column (14), 


Table 3 is set up so as to 
The start- 


ing conditions are Wy' = 1; Ww" = gio! = gw = 0. 
In Section A, the value 282.8 [= (1/2)k?Q?X?m,,_ (1/4) | 
is introduced in column (10), row z = 10, to conform 


with Eq. (66c). Similarly, the value 280.0 [ = (1/2) 
RQX?a,,—(1/4)] is introduced in Section B, column (14), 


row 7 = 10, to conform to Eq. (67c). 


y' = +0.00909378: 

Z,' = 1-0. $00061124 
yl = +0.457712776 
Z,1! = +0.919709104 


Step 5. Compute the Induced Moments.—The mo- 
ments at any station, 7, may be readily computed by 
the first two equations of Eqs. (65) and the 
previously evaluated. 


use of 
coefficients and parameters 
Table 4 shows this computation. 

Since the computations were rather lengthy, it was 
desirable to have a final check on the solutions. This 
was accomplished by means of a Check Table (not pre- 
sented), which is similar to that presented in reference 
3, with, however, the effect of the damping terms added. 
The solutions presented herein were checked, and the 
numerical accuracy was found to be very good. 


COMPARISON OF RESULTS AND CONCLUSIONS 


An approximate, as well as a more precise tabular, 
method for determining the bending moments induced 
in a rotor blade has been presented. Both methods 
include the effect of aerodynamic damping. In order 
to ascertain the influence of the aerodynamic damping, 
results were obtained for a typical helicopter blade in 
forward level flight and were compared with the re- 
sults obtained by neglecting this damping effect. The 
comparisons of results are shown in Fig. 2. The maxi- 
mum second harmonic station moments and their phase 
angles presented therein are the combined results of 
both the sine and cosine components of the second har- 
monic air loads. It is immediately evident that the 
damping effect is not a minor one, and, in the design of 
efficient blades, it should not be neglected. It is also 
noticed that good agreement exists between the mo- 
ments calculated by the approximate and more precise 
tabular methods presented in this paper. An appreci- 
able error does exist between the phase angles at the 
outboard end of the blade. However, at the important 
inboard stations of maximum moments, agreement is 
obtained. The disagreement near the tip results 
from approximating the deflections by a single mode 
shape. This approximation is reasonably good for the 
inboard stations but is in error at the outboard stations. 
However, the accuracy of the approximate solution could 
be improved by use of the more accurate method involv- 
ing carrying out Step 6, as indicated previously. 

It is therefore concluded that the amplification fac- 
tor method, including aerodynamic damping, may be 
used to quickly and accurately determine the maxi- 
When, however, greater 


mum rotor blade stresses. 
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accuracy is desired, one may resort to the tabular 


method. 
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Appendix 


SANDWICH PLATES WITH ANISOTROPIC CORE 


Although not published in reference 1, it may be use- 
ful to give here the necessary data for sandwich plates 


with anisotropic core. This affects only P, for Case 2. 
Denoting the moduli of rigidity for deformation by the 
transverse shears Q, and Q, by G,, and G,,, the differ- 
ential equation (81) for Case 2 transforms into 


(t + h)? O*w O7w 
Ges ( ° + € ‘) ee 
t Ox? Oy? 


O7-w O*w "Ww 
i= - fF, — + 3 Pe =(0 (91) 
Ox? Oy" Oxoy 
where 
e = G,,/Ga (92) 


In the same way as for isotropic core one obtains, in- 
stead of Eq. (30), 


P, = al(t + h)°G,,/t] (93) 
a is given in Table 3, where again 8 = L/b and y = 
Fist ws 
TABLE 3 
Case a 
A 1 
Bf, D 1 + ef? 
E, F 1 + (€8?/4) 
G.H [e(1 + €82)}! 
I (1 + €8?)/(1 + 6?) 
K e| —(y/2) + [(y?/4) + (i/e) + 8B?) 


All other formulas and values in Table 1 remain the 
same, except that in the expression for 7 now G, has to 
be replaced by G,,. Furthermore in Case I, assuming 
the plate to have infinite length, ky) = 1 + (1/8°). 











The Classification of One-Dimensional Flows 
and the General Shock Problem of a 
Compressible, Viscous, Heat-Conducting 


Fluid’ 


GEOFFREY S. 


Harvard [ 


SUMMARY 


Recently a great deal of attention has been focused on the prob- 
lem of determining the structure of a shock wave in a perfect gas 
on the basis of the Navier-Stokes equations. In this paper, we 
shall denote by ‘‘perfect gas’’ a viscous, heat-conducting fluid 
that satisfies the equation of state p = 
are in general functions of temperature. 
consider the exact structure for a case where the Prandtl Number is 
His analysis was restricted to constant 


gRpT, whose specific heats 
Becker? was the first to 


neither zero nor infinite. 
specific heats, constant coefficients of heat conduction and vis- 
cosity, and constant Prandtl Number = 3/4. After him several 
authors*® ®-§ have attempted to solve the problem under 
more general conditions. Since the first approximate discussion 
by Taylor,'! other approximate methods have been developed.* !! 
In some cases, the Navier-Stokes equations were reduced to an 
ordinary differential equation of higher than the first order; in 
others the reduction of Becker to a first-order equation was used. 
The approximate methods resort directly to the Navier-Stokes 
equations. 

It was von Mises’ who reformulated Becker’s approach and 
thus gave not only a much simpler presentation of the shock 
problem, but also a method of reviewing the general charac- 
teristics of all possible one-dimensional steady motions. It is 
the purpose of this paper to show that, if the Navier-Stokes equa- 
tions for the flow of a perfect gas are taken as a basis, then the 
method of review given by von Mises enables a complete classifica- 
tion of the flows to be made. Moreover, a discussion is presented 
of the problem of the shock wave in more than one dimension, 
which shows that a fuller numerical treatment of von Mises’ equa- 
tions and not a new!! approximation is needed in this case. 


CLASSIFICATION FOR CONSTANT SPECIFIC HEATS 


| 2 pee THE SOLE ASSUMPTION that the specific heats 
of the perfect gas considered are constant, von 
Mises® reduces the Navier-Stokes equations of the 
motion to 


dt I 4P t - to =! 4P 


= - at. <5 << oO (1) 
dv 3Y t —_ oe 3dY 
where 
t = gRT/C,’; vy = u?/2C,? 
bob = (y-lw-vV2a+ 0); to = V 20 — Ww 
Cc = Cs Ci? 
Received March 14, 1951. 
* Paper done under Contract Ndori 76/XVI NR 043-046, 


with the Office of Naval Research. 
+ Teaching Fellow in Applied Mathematics. 


S. LUDFORD? 


‘nwersity 


and C;, Cy are certain energy and momentum con- 
stants. This equation may be discussed in a qualitative 
manner so as to give a classification of all possible flows 
without recourse to numerical integration. For the sake 
of clarity, we shall do this for the case P constant and 
equal to 3/4; this case not only has the unique property 
of possessing an integral* of the same form as Ber- 
the only known closed solution, first 
but also exhibits properties discussed 


noulli’s equation 
found by Becker 
by other authors. 
different constant is represented in the (v, ¢) plane by a 
the integral curves f 


The case of P variable or equal to a 


geometrical distortion of for 


P = 3/4. 


The integrals in the (v, ¢) plane of the ordinary dif- 
ferential, Eq. (1), have singular points 1 and 2 at the 
two intersections of the parabolas 


(2a) 
(2b) 


The coordinates of these points satisfy the Rankine- 
Hugoniot conditions for a shock discontinuity. The 
two pairs of slopes of the integral curves at these points 
are given in von Mises’ paper, and the general shape of 
the curves is as shown in Fig. 1. The normal shock solu- 
tion is 1. In addition to these curves through the 
singular points, integral curves of all possible types are 
shown. We note that IV cannot intersect the sonic 
line since its slope at 1 is less than that of the latter. 
It seems unlikely that IX will either, except for J, 
nearly equal to one. Further, since for large v, ¢ Eq. (1) 
takes the form 


[¢ — (vy — 1)v]/(t + 22) } 
(P = 3/4) 


dt/dv = (4P/3y) } 


(3) 


we see that the integral curves tend to have asymptoti- 
cally either the slope —1 or —(y — 1)/y: thus XVII, 
XVI, XIII, IH, III, and IV all eventually cross the v- 


* Hereafter this integral will be referred to as Bernoulli’s equa- 


tion 
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axis. A further discussion of Eq. (3) shows the latter 
statement to be true whenever P is less than 37/2.* 

There is a second case to consider (see Fig. 2)—i.e., 
when the points | and 2 (hereafter called the Rankine- 
Hugoniot points) are imaginary: this is the case when 
c> y*/|2(y? — 1)], where c is a constant introduced by 
von Mises [see Eq. (1)]. In this case also, the curves 
have the asymptotic behavior of the last paragraph. 

Our interest is clearly limited to the first quadrant of 
the v, ¢ plane, since, on the v-axis, 7’ becomes zero, and 
on the ¢-axis, p becomes infinite. Moreover since Eq. 
(1) is derived from the two equations 


(4/3)(u/m)(dv/dx) = t —t (4a) 
[(y — 1)/gR](k/m)(dt/dx) = t — to (4b) 


and m is a positive constant, the direction in which x in- 
creases along an integral curve is determined uniquely 
at each point and is indicated in the figures by arrows. 
As we approach the R-H points, x > * @. 

* Here P is considered constant, but a similar criterion applies 


when it is variable. 
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In each case, it will be interesting to discuss the 
change in entropy along the flow. Here we must avoid 
some confusion in the literature.’ The “shock wave’”’ 
solution represented by the singular integral curve 
joining | and 2 in Fig. | is the only flow that can be 
maintained from x = —o tox +; all other flows 
must either be started or terminated, or both, at finite x, 
since otherwise v or /, or both, would be negative for part 
of the flow. In order to be started or terminated in a 
manner that satisfies the boundary conditions, some ex- 
ternal agency is required. Thus, the fact that the 
entropy may have a smaller value at the finish of such a 
flow than it has at the beginning is not an indication 
that the flow is physically impossible. The Second Law 
of Thermodynamics applies to a closed system. How- 
ever, the physical significance of these terminating 
solutions is open to conjecture. 

We have for any two points in the (v, ¢) plane on the 
same integral curve: 


S] = [gR/(y — 1) ]flog Tu?~*] = 
[gR/(y — 1)]flog w'?~ 7] (5) 





where [| | denotes the increment. Thus, all points 
(v, t) which have the same entropy as the point (v, fy) 
and can be reached from the latter along an integral 
curve, lie on 

S = ty ~ 9)? = ty,~-?” (6) 

It will be seen from Figs. 1 and 2 that there are three 
divisions into which the flows fall: 

(a) Those originating at x = —o© with uniform 
supersonic conditions. These are represented by the in- 
finity of curves leaving point 1 in Fig. 1. Points on 
these curves having the same entropy as | lieon S = S$). 
The slope of the latter at 1 is 


(dt/dv), = —[(y — 1)/2](t4/m) = 
—(y — 1)[(1/v2,) -—1] (7) 


since the point 1 lies on (2b). Hence, it touches (2a) at 
the point 1. 

(b) Those originating at x = — © with uniform sub- 
sonic conditions. These are represented by just the two 
curves XII and XIII leaving point 2 in Fig. 1. Points 
on these curves having the same entropy as 2 lie on S = 
5S». Clearly the latter touches (2a) at the point 2. 
Moreover, it passes above 1 since it is known that 
entropy increases through the shock. An immediate 
deduction is that the entropy has a maximum some- 
where within the shock; this is true whatever P may be 
(even variable). For P = 3/4, Morduchow and Libby’ 
have placed this maximum at the center of the shock. 

(c) Those originating from rest conditions at finite x. 
These are represented by curves starting at the /-axis 
and occur in both figures. Points on these curves having 
the same entropy as at the starting point lie on one or 
other of the axes. 

The results of applying the above analysis to Figs. 1 
and 2 are set out in Tables 1 and 2. We use there the 
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TABLE 1 
Flows Originating from Uniform Conditions at x = — @ 
(a) Supersonic atx = —« Reference Entropy 
I Ordinary shock-wave solution. Flow proceeds to uniform 2,3, 5,7,8,10 Entropy reaches a maximum for some finite x 
conditions at x = +o, Bernoulli’s equation satisfied for 
P = 3/4 
II Compression wave followed by expansion wave. Flow ter- 5, 7, 10 Entropy has a maximum in the flow 
minates at finite x. Partially subsonic. JT has a maximum 
III Similar to II except that flow is always supersonic Entropy has a maximum in the flow 
IV Expansion wave, terminating at finite x. 7 has a maximum Entropy has a maximum in the flow 
V_ Expansion wave, terminating at finite x. JT has a maximum Entropy eventually has value less than at 
x=- 
VI Expansion wave, terminating at finite x. Bernoulli’s equation 7 Entropy eventually has value less than at 
satisfied for P = 3/4 x = —o, (For P = 3/4 entropy decreases 
monotonically ) 
VIL Expansion wave terminating at finite x Entropy eventually has value less than at 
x =-a 
VIII Initially an expansion wave, then a compression wave; ter- ; Entropy eventually has value less than at 
minating at finite x x=-—« 
IX Compression wave terminating at finite x. 7 monotonically Entropy decreases monotonically 
decreasing 
X Compression wave terminating at finite x. Wholly supersonic. Entropy has maximum in flow 
T has a maximum 
XI Compression wave terminating at finite x. Partly subsonic. 5, 7, 10 Entropy has maximum in flow 


T has a maximum 
(b) Subsonic atx = —«@ 
XII Compression wave terminating at finite x. 
decreasing 


XIII Expansion wave terminating at finite x. 


expressions “‘compression wave”’ and ‘‘expansion wave” 
to denote that p increases and decreases in the direction 
of the motion, respectively. The references give the 
papers where explicit numerical integration has been 
recorded. 


CONCLUSIONS 


A certain number of conclusions may be drawn from 
a survey of the figures and tables. 
(1) If P = 3/4 and Bernoulli’s equation is satisfied 


with subsonic conditions at x = — ©, then the flow is 
uniform. 

(2) If P = 3/4 and Bernoulli's equation is satisfied 
with sonic conditions at x = — ~, then the flow is either 


uniform or follows VI. 

(3) The “shock wave” is the only nonterminating 
flow. 

(4) We see the reason for the instability in the 
numerical calculations of the shock-wave flow found by 
Reissner and Meyerhoff.© This was explained by 
Meyerhoff on the basis of a paper by Morduchow 
and Libby.’ It is clear from Fig. 1 that a small error in 
determining a point on the curve I between 1 and 2 will 
result in a divergence of the numerical solution along a 
curve of type IT, III, X,or XI. The numerical sensitivity 
of curve I to small error is thus noted for arbitrary P. 
Moreover, it is clear that the linearized theory carried 
out by Reissner and Meyerhoff*® suffers the same draw- 


T monotonically 


T has a maximum 


Entropy decreases monotonically 


Entropy has a maximum in the flow 


back. We should, in fact, expect any exact or approxi 
mate integration starting from the point | to give a fair 
representation only over the first part of the curve I. 
However, a backward integration from the point 2 
would not exhibit this instability, since there are only 
two curves passing through 2 in each of the singular 
directions (I, XIV; XII, XIII). 


(5) We see from Eq. (31) of von Mises’ paper that, for 
P < 3/4, the shock-wave integral lies everywhere 
nearer to, and, for P > 3/4, everywhere further from, 
the parabola ¢ = ¢., (this parabola is fixed for whatever 
P is, as also is ¢ = 4))._ Hence, from Eq. (4a) we see that 
for the same m 


P < 3/4: u(dv/dx) is everywhere smaller in 
magnitude 
P > 3/4; u(dv/dx) is everywhere larger in 


magnitude 
than for P = 3/4 (for the same v). However, since the 
change in uw consequent to the change in P is not pre- 
scribed, we cannot in general say anything about the 
effect on the velocity gradient.* ! If uw is unaltered, we 
may say the velocity gradient is smaller (or thickness 
, 


larger) for P < 3/4. 


It may again be stressed that although the above 
classification was carried out for P = 3/4 (u, k variable), 


the conclusions are easily generalized. We merely have 
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TABLE 2 


ONE-DIMENSIONAL FLOWS 833 


» 


Flows Originating from Finite x from Rest Conditions 


(c) Real R.H. points 


XI\ Expansion beginning at finite x and ending in uniform conditions at x 
noulli’s equation satisfied for P = 3/4. Wholly subsonic 


X\ Expaasion beginning at finite x and terminating at finite x as slight compression wave 


7 monotonically decreasing 


XVI Expansion beginning at finite x and terminating at finite 


XVII Expansion beginning at finite x and terminating at finite a 


(d) Imaginary R.H. points 


XVIII Similar to X\ 


XIX Expansion beginning at finite x and terminating at finite x 


XX Expansion beginning at finite x and terminating at finite x 


satisfied for ? = 3/4 
XXI_ Similar to XIX 
XXII Similar to XVI 
XXIII Similar to XVII 


a distortion of the integral curves in Figs. | and 2 with 
the A-curves preserved 


VARIABLE SPECIFIC HEATS 


This case has been considered by many au 


thors.* 46% 1° The analysis of von Mises may be 
modified by substituting in Eq. (3) of his paper 
oo Si is gR , 
K(T) C, dT for 1 (S) 
7 0 ¥Y l 


Thus, instead of Eq. (1) we have 
dt/dv = K }fe(t) — to*|/(t — ta)§ (9) 
where 
e(t) (1/C,), E(C7t/gR); K = 4ugR/3k 
to* T~-~ Vate 


Once the graph of £(7) versus 7 has been obtained, the 
graph of e(t) versus ¢ can be determined for a specific C; 
by a suitable stretching of ordinates and abscissas. The 
graph of A(t) versus ¢ can be determined likewise. 

The two singular points of the first-order equation, 
Eq. (9), are given by the intersection of the curves 


t—t. =0 (10a) 
e(t) — &* = 0 (10b) 


and the problem of determining the shock solution is 
that of determining the singular curve joining these two 
points. 

The advantage of this approach is that not only is the 
ordinary differential Eq. (9) an equation of the first 
order, but also its integrals admit the simple geometrical 
representation. In fact, in the preceding classification, 
we have an excellent guide to its behavior. In the other 
discussions in the literature, both these advantages are 
lost. In particular, we may expect that “backward” 
integration will prevent instability of the numerical cal- 


culations. 


Entropy 


éntropy monotonically increasing 


4 


maximum in the flow 


tntropy 


maximum in the flow 


x. TJ initially decreasing Entropy has 


maximum in the flow 


T initially increasing Entropy has 


éntropy has a maximum in the flow 


maximum in the flow 


éntropy has 


maximum in the flow 


Jernoulli’s equation Entropy has ; 


tntropy has a maximum in the flow 


entropy has a maximum in the flow 


Entropy has a maximum in the flow 


FURTHER ASPECTS OF THE THEORY 


In a recent paper,'' an approximate method is applied 
to the determination of the shock-wave solution for con 
stant gas coefficients. In essence, the authors introduce 
an approximation into the fundamental Navier-Stokes 
equations. Thus, instead of the energy equation, a 
polytropic pressure-density relation 


p = Ap’ (11 


is substituted, the other equations being unchanged. 
With this approximation, the equations can be dealt 
with more simply, and it is found that » must be a cer- 
tain function of the initial Mach Number in order that 
the Rankine-Hugoniot conditions are satisfied at 
=a + oe, 

The authors find on the basis of this arbitrary simpli- 
fication of the basic equations that the resulting shock- 
wave solution is a close approximation to the exact 
numerical one of other authors. They therefore suggest 
that such a simplification in the case of the equations of 
multidimensional flow will ensure mathematical tracta- 
bility without loss of close approximation. 

There are several objections to this: 

(1) It is not evident either physically or mathe- 
matically what is the nature of the approximation. 
Moreover, it would be more surprising if there were not 
close approximation in the one-dimensional case. 

(2) Even after the simplification has been made, the 
resulting equations must still be solved, and it cannot be 
foreseen that this is essentially simple. 

(3) Since m depends on the initial Mach Number in 
the one-dimensional case, we must expect it to depend 
in an even more complex manner on the flow outside the 
shock transition region in the multidimensional case. 

Against this may be put a treatment of the problem 
developed by the author.'” By applying an analysis 
that is completely analogous to that used in obtaining 
the Prandtl boundary-layer equations, a new set of 
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“shock transition’ equations can be set up, integrals of 
which represent asymptotic solutions of the Navier- 
Stokes equations as the Reynolds Number tends to 
zero. In the same way that the difference between a 
solution of the boundary-layer equations and the corre- 
sponding solution of the Navier-Stokes equations is 
known to be of 0(R7' 
set of “shock transition’ equations is known to be of 


O(R-?). 


*), so the difference for this new 


Furthermore, these ‘‘shock transition’? equations are 
identical in form with the one-dimensional equations 
discussed by von Mises® and can be obtained for arbi- 
trary gas coefficients. Indeed, they are of the same 
form for unsteady motion. Thus, the whole multi- 
dimensional problem is reduced immediately to a corre- 
sponding one-dimensional problem, the type and extent 


of the approximation being well known. 


It may therefore be asserted that what is required in 
the treatment of the general problem is not a method of 
rendering the basic equations more simple, but a com- 
plete tabulation of shock solutions in the one-dimen- 
sional case under change of arbitrary constants. 
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Horizontal Tail Loads in Abrupt 
Pull-Ups from Level Flight 


JAMES L. 


The Glenn L. . 


ABSTRACT 


A method is presented for the calculation of horizontal tail 
loads in abrupt pull-ups from level flight. The method is also 
applicable to pushover maneuvers. The basis of the procedure 
is the rational specification of the variation of load factor with 
time in the maneuver. The elevator deflection required to pro 
duce this airplane response is then calculated for the determin- 
ation of the maneuvering tail load 

The airplane response is shown to be a function of the maximum 
elevator deflection rate in the maneuver. This becomes the sole 
arbitrary variable in the method and has been tentatively cor- 
related against the airplane limit load factor and flight speed. 

The calculations are appreciably shorter than those entailed in 
existing methods, and, since only one arbitrary variable is pres- 
ent, it is believed that calculations for various airplanes can be 


carried out on a more consistent basis. 


INTRODUCTION 


i ip DETERMINATION OF THE HORIZONTAL TAIL loads 
in maneuvering flight is of significant importance 
because critical tail loads generally result from maneu- 
The total horizontal tail load 
One part 


vering flight conditions. 
may be considered as the sum of two parts. 
is the balancing load required to trim the airplane in the 
level flight attitude. The second part is the maneuver- 
ing tail load that is the incremental tail load developed 
as the airplane load factor is altered in an abrupt ma- 
neuver. The flight maneuver to be considered in an 
analysis to determine the critical horizontal tail load 
involves motion in the plane of symmetry only. The 
maneuver is of short duration so that pitching acceler- 
ations requiring large tail loads are developed. 

The basic problem in finding the horizontal tail load 
is the determination of the elevator position and tail 
angle of attack throughout the maneuver. Once these 
are known, the tail load is easily computed. 

Kelley and Missall' have developed a method for the 
computation of the maneuvering tail loads which has 
been accepted as the standard method of computing 
these loads by the Air Forces. Their procedure is to 
determine the dynamic response of the airplane to 
arbitrary elevator motions and to alter the arbitrary 
elevator motion until a satisfactory airplane response, 
as indicated by reaching a specified load factor within 
a particular time interval, is obtained. The solution 
to this problem is not unique in that many different 
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elevator motions will attain this end. The computa 
tions involved are also comparatively lengthy. 

The method proposed in the present report attacks 
the problem in an inverse manner. The airplane re 
sponse is specified in a rational manner and the elevator 
motion required to obtain this response is determined. 
The solution thus obtained is unique, and the use of 
but one arbitrary variable, the maximum rate of ele 


vator deflection in the maneuver, is required. 


SYMBOLS 


L = lift force, lbs 

M = pitching moment, ft.lbs. 

| = airplane pitching moment of inertia, slug ft.” 
Z = resultant force normal to X axis, taken posi 


tive downward, lbs 
Ly = total horizontal tail load, lbs. 
= tail load required to produce an increment of 


load factor, Ibs 


La, = tail load due to horizontal tail angle of attack, 
Ibs. 

Li. = tail load due to elevator deflection, lbs 

S = wing area, sqft. 

V = airplane true air speed, ft./sec 

V és = airplane true stall speed at the limit load factor, 
based on static CL,,,,, ft./sec. 

a, b, c, d, e = constants in the load factor-time equation 

c = wing mean aerodynamic chord, ft. 

1 = tail incidence relative to wing chord line, rad 

l = tail length, ft 

m = airplane mass, slugs 

g = acceleration due to gravity, ft./sec.* 

n = airplane load factor in g units, normal to Y axis 

Cmax = airplane limit load factor, g units 

q = pitching velocity of airplane X axis, rad./sec. 

q = free-stream dynamic pressure, lbs. /sq.ft 

t = time, sec. 

tma = time required to reach the limit load factor, sec. 

Ta = vertical translational velocity along Z axis, 


positive when downward, an increment from 


level flight condition 


a = angle of attack, rad 

Aé, = change in elevator deflection from trimmed 
position, rad 

€ = downwash angle, rad 

7 = pitch angle between flight path and horizontal, 
rad 

0 = pitch angle between wing chord line and the 


horizontal, rad 
n = tail efficiency, ratio of dynamic pressure at tail 


to that in the free stream 


Cre = lift curve slope dC /da, per rad. 
Cua = pitching moment curve slope, dC,,/da, per rad. 
+ ee = maximum lift coefficient of airplane 
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Fic. 1. Nomenclature diagram. 

de/da = rate of change of downwash with wing angle of 
attack 

Cm, = rate of change of pitching moment coefficient 
per radian elevator deflection 

Co = rate of change of pitching moment coefficient 
with angular rate of change of angle of 
attack, per rad. per sec. 

Cg = rate of change of pitching moment coefficient 
with angular rate of change of pitch angle 
0, per rad. per sec. 

dei /ad, = ratio of elevator effectiveness Ch, to stabilizer 
effectiveness Cie, 

Cus, = rate of change of horizontal tail lift coefficient 
with elevator deflection, per rad. 

Dots above a quantity indicate successive derivatives with 


respect to time—e.g., 6, = dé./dt; 6 = d?0/dt? 


Subscripts 
w = wing 
t = horizontal tail 
e = elevator 
cm = complete model 
LF = level flight condition 


DEVELOPMENT OF METHOD 


The nomenclature diagram defining the axis system 
used in the present analysis is shown in Fig. 1. The 
axes are fixed to the airplane and rotate with it. The 
translational velocity along the Z axis, w, is zero in level 


flight. In the maneuver, 


w= VAa 


where Aa is an increment from the level flight angle of 
attack. 


It is customary in the computation of the maneuver- 
ing tail loads to assume that the air speed is unchanged 
over the short time interval that is required to perform 
the pull-up. This time interval is on the order of 1 to 
2 sec. With this simplification, the problem is reduced 
to one involving two equations of motion. 


Kelley and Missall,! in their investigation, have 
shown that certain of the terms in the remaining mo- 
ment and normal force equations of motion may be 
neglected with small loss of accuracy in the tail-load 
calculation. Eqs. (1) and (2) formed the basis of their 
method and also are the foundation of the procedure of 
this paper. 
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m(w — V6) = w(dZ/dw) (1) 
Ty, = w(0M/dw) + w(0M/dw) + 6(0M/d6) + 
A6,(0'M/06,) (2) 
Eq. (1) may be solved for 6 to give 
6 = w/V — [w(dz/dw) |/mV (3) 


The resultant force normal to the X axis of the air- 
plane has its positive direction downward along the Z 
axis. This is essentially opposite the lift, and we may 
write 


Z= -L 
w(0Z/Ow) ~ —Aal(OL/Oa) =~ —AaChacygeS (4) 


Now 
w= Va 
and 
w= Va 
Therefore, 
6 = & + [AaCLacygoS(g)|/[mV(g)] (5) 


The level flight lift coefficient 
(Cipp = mg/(goS) (6) 
Also 
Cipyp = Clacy Gir (6a) 
Substitution of Eqs. (6) and (Ga) in Eq. (5) yields 
6 = & + (g/V) (Aa/azr) (7) 


It is convenient to work in terms of the load factor, 
which we shall define as 


a= L/W 


It is easily seen that this relation may be written, as 


n = A/ Qype (S) 
From this, 
(n — 1) = Aa/ay (9) 
and 
n= a&/apr (10) 


Eq. (7) may now be written as 
6 = apn + (g/V) (n — 1) (11) 
6 = appt + (g/V) 0 (11a) 


Utilizing the above equations, Eq. (2) may be reduced 


to coefficient form and stated as follows: 


» fey e g lis ce 4 c | 
nN , n _ A < : 
LF quSé V qoSé Apr\em M5 +t 


: £ an — 
ii 1) | — aux M . V c “| = A6.CMs, (12) 
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HORIZONTAL 


Eq. (12) expresses the dynamics of the pull-up from 
level flight in terms of the load factor and its time de- 
rivatives. Therefore, a specification of the time vari- 
ation of the load factor will enable the determination of 
the elevator position throughout the maneuver. 


DETERMINATION OF LOAD FACTOR vs. TIME CURVE 


During the maneuver, five boundary conditions exist: 

(1) m = land y = Oat¢t = 0: In this analysis, t = 0 
at the instant the elevator motion initiating the ma- 
neuver is begun. 

(2) m = Oatt = 0; 
steady flight at the beginning of the maneuver, there 
is no time rate of change of the load factor. 

(3) 2 = Oatt = 0; the pitching moment equilibrium 
existing provides no source of pitching acceleration. 

(4) m = Mmar. Att = tmar.3 by definition, since tmaz 
is the time required to reach the maximum load factor. 

(5) 2 = Oat tmar.; the derivative of a function is zero 


since the airplane is trimmed for 


atits maximum. 

These five boundary conditions enable the deter- 
mination of five constants for a mathematical repre- 
sentation of the » ~ ¢ variation in the pull-up. Studies 


N( ape(Lyy,/goSé)| + #[(g/V) Lyy/qoSé) — ape(Cr = is Cy;) | v h{—aprCu, 
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indicated that a polynomial expression provides a satis- 
factory relationship, as will be shown later in the 


paper. 
n=a-+ bt + ct? + dt? + et 
From condition (1), a = 1; from condition (2), b = 0; 
from condition (3) ¢ = 0; from conditions (4) and (5), 
d = [4A(Mmaz. Eh, So 
and 


é = [—3(mar ohh 


Thus, the resulting expression for the variation of 
load factor in the pull-up is 


i= ] + [ (mar -_ RD hice *1 x 
[4t® — (3f*/tuer.) | (13) 


SOLUTION OF DyNAmIcC PULL-UP EQUATION 
If Eq. (12) is differentiated with respect to time, the 
resulting expression provides a method of determining 
the rate of elevator deflection throughout the maneuver. 


(g 1) Carg | = 5.Cur,, (14) 


Using Eq. (13) and its time derivatives with Eq. (14), it is possible to express the maximum rate of elevator de- 
flection in the pull-up maneuver in terms of the airplane stability and inertia parameters, the elevator control 


power, the maximum load factor, and the time required to attain this load factor. 


a 
_—_ a E : 


employed, the following expression results: 
— 12( maz. = i) J I 


emar - ~ 
é Mse(tmaz. )8 


It should be noted that the maximum elevator de- 
flection rate given by Eq. (15) is 6, for the time at which 
Nmar. iS attained. This is a consequence of the vari- 
ation of load factor with time given by Eq. (13), which 
requires é, to be a maximum at faz. 

When a value is specified for é, mar. the maximum 
rate of elevator deflection during the maneuver, it is 
possible to solve for the complete time history of 
maneuver by means of the methods of this paper. The 
time to attain the maximum load factor may be found 
by solving Eq. (15) for tinu:.. 

The terms (Cy; + Cu.) in Eq. (15) represent the 
damping in pitch of the airplane. It 1s customary to 
calculate Cu, by increasing the direct contribution of 
the horizontal tail by 25 per cent. This allows for the 
contribution of the other components of the airplane 


to the damping in pitch. Cy, results from the lag 


When this procedure has been 


: . { oe 
qoSE — appr(C Mg + C “|| (15) 


in downwash at the tail surfaces. (Cy, + Cy.) are 
a 


then evaluated as 


(Cuy + Cur,) = Cmi(l, V) [(1.25 Vv n) + (de da) | 
(16) 


CALCULATION OF THE MANEUVERING TAIL LOAD 


It is now possible to determine the tail angle of attack 
and elevator deflection throughout the pull-up ma- 
neuver and thus to find the maneuvering tail load. The 
change in tail angle of attack is given by the formula: 

Aa, = Aay — Ae + O(1,/ Vn) (17) 

The downwash angle « should include the lag in 

downwash which results from the finite time interval 


required for the air to flow from the wing to the tail. 


This lag angle is calculated as 
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€lag = —azrn(l,/V) (de/da) (18) 


The total downwash increment is the sum of the 
downwash lag angle and the downwash change that 
would exist due to wing angle of attack under static 


conditions. The total downwash angle change is 


Ae = [(m — 1)azp(de/da) — [azen(l,/V) (de/da)| (19) 


We then for the 


change in tail angle in the maneuver: 


have the following expression 
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Aa; = (nN — l)azr[1 one (de da) | + aprn(l, V) x 
[((1/-Wn) + (de/da)| + [(gl,)/(V2Wn)] (n — 1) (20) 

The increment in elevator deflection required in the 
pull-up is influenced by three factors: (1) the static 
longitudinal stability, (2) the damping in pitch, and 
(3) the pitching inertia and acceleration. The com- 
bination of Eqs. (12) and (16) yields the following ex- 
pression for the increment in elevator deflection from 
the level flight value. 


P= appl yy g | in QLrF i ] 2D de 
Ab, = %| = —-|+n|/—— — : = + + 
CM5,qoS€ V CM5.q0S€ (da,/db.) V\ V9 da 


The maneuvering tail load is calculated from the 
equation 


Doce. = gon S ,[ Aa,Chay + A6,CL,;, | (22) 


DISCUSSION 


Fig. 2 shows the variation of elevator deflection and 
load factor with time in a pull-up performed in flight 
with the JRM-1 airplane. 
calculated with the polynomial load factor versus time 


For comparison, curves 
relation are shown. Inspection of Fig. 2 shows that 
the recommended polynomial is a good representation 


of the flight curve. It may also be noted that the agree- 
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ment between the actual and computed elevator de- 
flections is good. In both cases, the elevator deflec 
tion rate is more rapid in the checking of the maneuver 
than in the entry. This is consistent with data on 
other aircraft which have been reviewed by the author. 

It is evident that the tail angle of attack is relatively 
near its peak at the limit load factor, since the wing 
angle of attack is at its maximum at this point. The 
reversed downward tendency in the elevator motion 
near the limit load factor adds to the upward tail load 
due to tail angle of attack. Therefore, the important 
problem in the determination of the maneuvering tail 
load is the prediction of the maximum downward ele- 
vator deflection at the limit load factor. The down- 
ward elevator tendency is required to produce the nega- 
tive pitching acceleration needed to check the maneuver. 
It is apparent, then, that the pitching acceleration 
existing at the limit load factor is a vital factor in the 
determination of the critical tail load in the maneuver. 
The rapidity of the maneuver dictates this pitching 
acceleration to a large extent. When Eq. (13) is used 
to compute the load factor-time curve, the pitching ac- 
celeration is uniquely determined throughout the ma- 
neuver. Therefore, it is of primary importance to deter- 
mine the time interval required to attain the limit load 
factor to be used in Eq. (13), since this time interval 
fixes the pitching accelerations for a particular set of 
airplane flight conditions. 


MAXIMUM ELEVATOR DEFLECTION RATE AND ITS 
EFFECT ON PULL-UP TIME 
It may be seen from Eq. (15) that the time required 
to perform the pull-up is a function of the maximum 
rate of elevator deflection. If a rational specification 
of this parameter can be determined, then the calcu 
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lation of horizontal tail loads for all airplanes may be 
accomplished on a consistent basis. In the specifica- 
tion of the maximum elevator deflection rate, it must 
be borne in mind that the horizontal tail strength should 
be compatible with the wing strength. There is little 
point in designing the empennage for air loads, which, 
if applied, would result in failure of the wing. In line 
with this reasoning, the author believes that the maxi 
mum elevator deflection rate should be a function of 
the limit load factor of the airplane. The importance of 
this may be seen by considering the airplane of limit 
load factor of unity. Any elevator deflection rate 
that causes a positive change in the angle of attack 
before an air-speed decrement occurs will result in the 
limit load factor exceeding one. Therefore, since the 
equivalent air speed is assumed to be constant through- 
out the maneuver, the allowable elevator deflection rate 
is zero. 

It is also obvious that the air speed should have an 
influence on the allowable elevator deflection rate. 
The elevator hinge moment and increment in airplane 
pitching moment vary as the square of the air speed. 
It is reasonable to assume that the time required by the 
pilot to exert a given stick force is independent of the 
air speed. Therefore, the elevator deflection requiring 
a fixed pilot effort will vary inversely as the air speed 
squared. Consequently, the elevator deflection rate 
will also be an inverse function of the square of the air 
This factor may be nondimensionalized by 
\’)*, where 


speed. 
making the elevator rate a function of (V,, 
l’,, is the stalling speed at the limit load factor for the 
appropriate configuration and altitude. 

With the above considerations in mind, it appears 
reasonable to specify the maximum elevator deflection 
rate in the form of 

5 = Rta, — 1)*(V,,/V)* 


emar 


where & and a are yet to be determined constants. The 
studies that have been performed thus far indicate the 


following equation is suitable. 


Dense. = SO Ulnar, — 1)* (V,, 


COMPARISON OF PULL-UP TIMES 


Kelley and Missall' have proposed the following re 
lation for the maximum acceptable time required to 
develop the limit load factor for use in structural calcu- 


lations. 


bitan: 2 CH ions 1) ]9-357 (24) 


A comparison between /,,,;, computed from this for 
mula and fas 
(23) in Eq. (15) has been made for three Martin air 
This 


calculated using the elevator rate of Eq. 
planes at |’, and at the airplane limit dive speed. 
comparison is shown in Fig. 3. 

It may be seen from examination of Fig. 8 that the 
pull-up time computed by the method of this paper at 


l’,, compares closely with that determined from Eq 
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(24) for airplanes (1) and (3), being slightly more rapid 
In the case of airplane (2), a longer time 
How- 


in both cases. 
is indicated at V,, with the method of this paper. 
ever, it is interesting to note that the pitching moment 
of inertia of airplane (2) is approximately 16 times that 
of airplane (1), while the product of the tail area times 
the tail length is but 1.7 times that of airplane (1 
Therefore, it would be expected that much larger pitch 
ing accelerations could be developed by elevator motion 
in airplane (1), and consequently, a much more rapid 
maneuver could be performed with this aircraft 

From Eq. (24), the pull-up time is not a function of 
the air speed in the maneuver. The present study indi 
cates that the pull-up time should become more rapid 
as the air speed is increased. At the airplane dive 
speeds, the pull-up time for all three subject airplanes 
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is well below that given by Eq. (24). This is a conse- 
quence of the decrease in az; resulting from higher air 
speeds. Because of smaller values of a,x, a given 
pitching acceleration will cause more rapid load factor 
changes at higher speeds. It appears that some ac- 
count should be taken of this effect in the computation 
of the maneuvering horizontal tail loads. 


COMPARISON OF CALCULATED TAIL LOADS 


In Fig. 4, « comparison between the results computed 
by Kelley and Missall' and by the method of the pres- 
ent report is presented. The C-54 airplane! was used 
for the comparison. The maximum maneuvering tail 
loads computed by the two methods agree well, being 
14,070 Ibs. by Kelley and Missall and. 14,980 Ibs. by 
the method of the present paper. This comparison 
appears satisfactory. With both methods, the maxi- 
mum negative pitching acceleration occurred at the 
limit load factor, with the slightly higher value of the 
present report accounting for the higher tail load. 

As the tail load computed by the proposed method is 
still increasing as the maximum load factor is reached, 
the problem arises as to whether the recovery from the 
pull-up might not provide even higher loads. In order 
to investigate this problem, the pitching acceleration 
existing at %mar., Which is the largest negative pitching 
acceleration during the maneuver, was assumed to re- 
main constant throughout the recovery. The total 
tail load in the recovery, calculated under this assump- 
tion, was always less than that at the maximum load 
factor. Study of flight-test data also indicates that 
the pitching acceleration in the recovery is markedly 
less than that at 2na7,, so it may be concluded that an 
analysis of the recovery from the abrupt pull-up is not 
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necessary in the determination of the critical tail load. 
This agrees with the results of Kelley and Missall as 
shown in Fig. 4. 

The proposed method may be equally well applied 
to push-over maneuvers to negative load factor by 
simply using a negative value for -_ in Eq. (28). 


CONCLUDING REMARKS 


The present paper has given another approach to the 
controversial subject of maneuvering horizontal tail 
loads in an effort to provide a rational solution to the 
problem and a solution that eliminates a considerable 
amount of the tedious labor involved in the calcula- 
tions. The specification of a rational load factor-time 
relation avoids the necessity of solving the differential 
equations of motion as such and transforms the prob- 
lem to the determination of the elevator motion re- 
quired to provide the specified airplane response. The 
specified airplane response is, in turn, a function of the 
value of the maximum elevator deflection rate in the 
maneuver. The maximum elevator deflection rate 
then becomes the sole arbitrary variable in the method. 
The uniformity of tail-loads calculation permitted by 
this fact is believed to be an important advantage of 
the proposed procedure. The recommended equa- 
tion defining the maximum deflection rate of the ele- 
vator has two constants whose tentative values are 
subject to verification by future studies and usage of the 
proposed method. 
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Appendix 


A sample calculation will be carried out as an example 


of the procedure of the present paper. To facilitate com- 


parison with the results of Kelley and Missall, the example airplane will be that used in their report,! p. 154. 


Model C-54 


cg. = 43% mac. Alt. = 15,000ft. 5S, = 1,457 sqft. 1, = 48.0 ft. 1/Vo = 1.262 
V, = 225m.p.h. G.W. = 62,000 Ibs. S, = 324.88 sq.ft. Ty = 560,730 slug ft.” 
Rese, 2H é = 13.64 ft. T.V.. = 0.785 Cimaz. = 1.5 
Cui, = —2.68 Cmacy = —90.0561 de/da = 0.40 Cia, = 4.10 CL, = 1.96 
Cuz, = —1.56n Clacy = 5.14 n = 0.834 
GS = V [2nmar.\ W/S)](Crmar.P) = 307 ft. per sec. = 210 m.p.h. (equivalent air speed) 
V.,/V = 0.933 
—_ = 35(2mar. — 1):* (V../V)? = 35.8°/sec. SPECIFIED 
: —12(" maz. = 2) j Tey E I yy \ 
6,, 4 ‘ tmax : oad LF C . + C, 5 j 15 ) 
naz CM5,t mas 3 ’ OLP goSé + x V qoSé apr MG mt) { (io 


oe led 


apy = (W/S)/(qoClacy) = 0.0637 
Tyw/(qoSé) = 0.217 
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Aa, = 
Aa, = 


Aé, = 





HORIZONTAL TAIL 


i] 3 
Cu, + Cug = Ci 7 [(1.25/Wn) + (de/da)] = —0.547 
Bonar = (11.52 /tmar.*) (0.0553 + 0.0516lmar.) REQUIRED 
| py . var, (rad./sec.) me (deg. /sec. 
0.6 4.56 261 
1 .o | .23 70.5 
1.4 0.535 30.7 
baw = Lis 
1 + (mer. — 1)/tmar.? [442 — (304/t mor.) 
1 + 2.62¢% — 1.4974 
7.8602 — 5.9603 
15.72 — 17.S8f° 
(g/V) (n — 1) + ayer 


= (0.0772(m — 1) + 0.0637n 


0.0772 + 0.06377 


LOADS 


(n — 1) farell — (de/da)|] + gl, V?V nt + 7 are(l./V) [A/V 9) + (de/da) | 


0.0479 (1 — 1) 4+ 0.010937 


; apr! yy g QuE t, 4125 de 
n . _ + i 35 = ad oe - + aa is — 1) — ies 
CM5,qoS€ Ve C5,qoSé da,/db, V\ V9 da 


—0.0107% — 0.0400% —O0.02S2(7 — 1) 
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V9 1.938 2.02 —().32 0.136 0.0666 —0.1040 —5.97 9, 5SO 
2 2.42 1.01 —(6§.82. —0.357 0.0784 —0.0094 —0.54 11,280 
32 2.20 8 —10.35 —0.660 0.0718 0.0675 3.87 10,320 


) 


CMacy 


CM;, 


$, 660 





2 sec. for the specified elevator rate of 35.8° 
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Use of the Local Mach Number in the 
Prandti-Glauert Method 


E. V. Laitone 
Associate Professor, University of California at Berkeley, Calif 
June 15, 1951 


THE SECOND-ORDER PERTURBATION EQUATION 
B* SUBSTITUTING THE PERTURBATION POTENTIAL @¢ defined by 
u= U + ¢z, v= dy (1) 
into the potential equation for two-dimensional isentropic flow 
and then neglecting all the third and higher order terms of @ and 
its derivatives, one obtains the second-order perturbation equa- 
tion 


2 —s 
[« — Ma?) — Ma* = o(1 +> Mt) + yy = 


2 
M «? py Pubew + O(¢*) (2) 
In deriving Eq. (2) one writes 
ae SOO : 
a? = aot +> [Ut — (U + or)? — oy] 


so that 
1 l ‘ Ee . 
= a 1+ (vy — 1)Mo*pz UV + O(g?) 


as long as 


bz | >> |b? + oy? | 
1>/ (7 — 1)Ma%(1/U) 


The expression for the characteristics of Eq. (2) show that as 


long as 


. a ‘ 
0<(~ Ma) >[ ot Fo(i + 2 Ma) + 
M * 4 3 
M co “py V (0) 


then Eq. (2) is of the elliptic type with imaginary characteristics 
The inequality given in Eq. (3) defines a free-stream Mach Num- 
ber \/.~ less than unity for which a local supersonic region with 
real characteristics first appears. 

Then, as long as the inequality expressed by Eq. (3) is satisfied, 
we may use the following iterative procedure for the elliptic-type 
equation expressing pure subsonic flow: Initially, we obtain 
the first-order approximation ¢; as a solution of 


(1 — M o*)\bire +t Piyy = YU (4) 


satisfying the given boundary conditions. 
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Then we may write the second approximation for an iterative 


solution of Eq. (2) either as 


2 1-1 
M~* Ale t J 9 Ma ‘Yon, Picr 1 Py, | 5 
(1 — M wo? )p»,, + o2, /= 


= ay eet. ‘ 
M ~? Vv 1 + 2 M w* }bi,b2r, + 01,02, (6 


In either case, if the given boundary conditions are satisfied, 
then ¢2 will be a second-order approximation that will also con- 
tain @, the original first-order solution. Although Eq. (6) 
should provide a better solution than Eq. (5), since its iterative 
process is essentially more rapidly converging, still Eq. (5) is the 
one generally used because its solution is simply the sum of the 
general solution of Eq. (4) added to the particular solution for 
the right-hand side of Eq. (5). 

However, for a numerical solution in a small local region, Eq 
(6) would be the logical iteration step for given fixed boundary 
conditions. In this case, the values for ¢; may even be assumed 
essentially constant in a small region. For example, one may 
write the first-order solution at any given fixed point on the sur- 


face of a slender body as 


o 


i = —(U/2) (C»/W1 — Ma®)' 

oi, = UO 
where C,,. is of the same order as @ and is the surface pressure 
coefficient for the incompressible flow solution of Eq. (4) when 
Mao = 0 


direction of the free-stream velocity LU’ 


The local surface slope 6 is measured relative to the 


Substituting Eq. (7) into Eq. (6) we obtain 
a 
(i -Me) + (14 = 


COMPRESSIBILITY CORRECTION FOR H1GH SUBSONIC SPEEDS 


Cos 
Mo :) ae? : x 
V1— M.? 


P2rr T Oryy = 21 «70d» , (8 


Now it is difficult to obtain a simple explicit expression for the 
surface pressure coefficient C, from Eq. (8) because of the ¢,, 


term and the fact that the boundary condition at the body sur 
face requires that 


— ... - <1¢ ; $s + Og» (9 
l T Pr l U 


However, two simplificatious may be logically assumed that will 
give a simple relation for Cp. (1) We apply Eq. (8) to high sub- 
sonic speeds so the right-hand side of Eq. (8) may be considered 
zero. This is a valid assumption because, at high subsonic 
speeds, 8 << 1, sothat 

(10 


(2) We assume that Eq. (9) may still be approximated to the 
second order by 

o2, = U6 (11 

It is easily shown that this is also a valid assumption for 8 << 1 

The second-order expression for the pressure coefficient is given 


c= -2(%) - (BY - (SY 


However, for a slender body, especially at high Mach Numbers, 


by 


only the first term is important in two-dimensional flow. There 
fore, since C,, is a constant for the given fixed point on the body 


surface, we obtain a solution of the simplified Eq. (8) as 
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Cp = —22,/U = Cp / Bx (12) 


where 


Sis 
2 = (1 — Mw?) + (1 = vet) Me? 
2 V1— Ma? 


If we reduce Eq. (12) so as to eliminate higher order terms of 
Cy, We obtain to the second order 
Ces 2M ao? + yMoao4 M 4 
C, = - —— Cy? (13) 
‘i — if.4 4(1 — M.?)? 

Eq. (13) is a new compressibility correction for two dimensional 
isentropic flow that should provide a satisfactory relation for 
thin bodies at high subsonic speeds. At low subsonic speeds 
Eq. (13) reduces to the Prandtl-Glauert relation, while at the 
higher subsonic speeds it predicts a slightly more rapid rise in the 
magnitude of the negative surface pressure coefficient than does 
the Karman-Tsien relation. The Karman-Tsien relation has an 
asymptotic value of 7 for infinite negative values of C, given by 

M.? = 1 (C,,/2)? 


which is practically unity for slender bodies. However Eq. (12 
has an asymptotic value of AJ. considerably less than unity and 
is in agreement with Eq. (3), the exact second-order expression 
predicting the transition from subsonic to supersonic flow in any 


given region. Therefore Eq. (2) is an upper bound. 


APPLICATION TO SUPERSONIC FLOW 


As a reasonable check on the assumptions made in obtaining 
Eq. (13) one may apply the same procedure to a pure supersonic 
flow where the known second-order solution is given by the 
Busemann relation. For pure supersonic flow the inequality is 
the same as that expressed in Eq. (3) if one writes 

(Mao? —1)>0 
The characteristics are now real, and Eq. (2) is always of the 
hyperbolic type. Then Eq. (7) is replaced by the well-known 
Ackeret solution 
i, = U6/V Ma? — 1 (14) 
where ¢;, is of course still defined by U@. Following the same 
approximations, Eq. (12) is replaced by 


C. = —2(g2,/U) = 20 By, (15) 


where 


y¥—1 26 
Br? = (Moo? — 1) (13 Mot) Ma’ 
2 V Mo? 1 


Then expanding to the second order in 6, one obtains 


2 2Ma* + yMa! — Ma! 
Cp = 64 . @? (16) 


The first term of Eq. (16) is, of course, the same first-order 
Ackeret solution. The second term containing 6? has the same 
denominator as the Busemann solution. The numerator is 
of a slightly different form, but as .1/. approaches unity the 
numerator approaches (1 + y), which is the form approached 
by the numerator of the Busemann relation. This is the result 
that should be expected, since the relation, corresponding to Eq 


10), 

(0/V Mao? — 1)d2,, >> O¢2,, (17) 
must be satisfied in deriving Eq. (15) from Eq. (6);  conse- 
quently, Eq. (16) is only applicable as 1/. approaches unity 

Use OF THE LocaL MAacuH NUMBER 


It is interesting to note that the local Mach Number 1/7 is 


always given in isentropic flow by 


2 VM «w? 
M,? = (1. + \(: + 
y-1 2 


so, if one defines 
Cp = —2(¢2/U) + 0(¢?) 


and expands Eq. (18) to the first-order, one obtains 


to 


M.«. :) Wao? — bi, + O(¢?) (19) 
U 

which is identical to replacing the expressions given in Eqs. (12) 

and (15) by 


Bp = V1 — M,? I 


(20) 
Bu= VMi? —1 § 


respectively. This is equivalent to using the local, rather than 
the free-stream, Mach Number in the ordinary Prandtl-Glauert 
correction. If Eq. (2) has the same properties as an ordinary 
linear equation, then the Prandtl-Glauert rule gives the lower 
bound, while Eq. (12) provides an upper bound for the increase 


of the peak Cp values with J/ &. 


Remarks on ‘The Decay of Isotropic 
Temperature Fluctuations in an Isotropic 
Turbulence”’ 


M. Z. E. Krzywoblocki 

Professor of Gasdynamics and Theoretical Aerodynamics, Secretary 
Panel on Fluid and Solid Mechanics, University of Illinois, 
Urbana, III. 

July 93, 1951 


§ pss IS A MISPRINT in the paper (see Stanley Corrsin, 
JOURNAL OF THE AERONAUTICAL SCIENCES, Vol. 18, No. 6, 
pp. 417-423, June, 1951). 
paper, the parameter c, should be used 


Instead of cp in y throughout the 
For an incompressible 
fluid this is true in both approaches: (a) the more elementary 
approach based on the mechanics of continua and (b) the 
approach based on the kinetic theory of gases (see S. Goldstein, 
Modern Developments in Fluid Dynamics, Vol. II, p. 606, and 
Chapman and Cowling, Mathematical Theory of Non-Uniform 


Gases). 


The Energy Equation for Two Kinds of 
‘*Incompressible Flow”’ 


S. Corrsin and L. S. G. Kovasznay 

Department of Aeronautics, The Johns Hopkins University 
Baltimore, Md. 

September 4, 1951 


ROFESSOR KRZYWOBLOCKE has pointed up the fact that refer- 
ence 2 was too vague in not designating the type of ‘incom 
pressible fluid’’ under discussion; therefore, some clarification is 
in order 
Two common types of ‘incompressible’ fluid are (a) the per- 
fect liquid (p = constant) and (b) very low Mach Number flow 
of a perfect gas 
The general energy equation can be written 
DT , OMe do 4, o7 
4+ p 


i = k + & 1 
Dt Ox; ox, | ox, | 


pe 
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where ® is rate of generation of heat by viscous dissipation and 
the other symbols have their conventional meanings. For a 
perfect liquid Ou;/Ox; = O, and, if the flow is slow enough to per- 
mit neglecting @ (~ velocity squared), 
pc(DT/Dt) = (0/Oxx) }R(OT Ox,){ (2) 

For low Mach Number gas flow things are more complicated. 
In this case, Eq. (1) can be transformed to 
DT Dp 0 S, oT | 
Dt Dt axe | Ox; f 
Now the question of whether pc,(D7/Dt) or pc,(DT/Dt) will 


be the correct approximation for the left side when 


&<X (0/dx;) | R(OT/dx;)} 


pCp + > (3) 


on the right side, depends upon the relative magnitudes of 
p(du;/Ox;) and Dp/Dt. 

As pointed out recently by Goldstein,‘ the correct choice prob- 
ably depends upon the particular type of flow. It appears that 
the correctness of neglecting Dp/Dt can be shown for a very gen- 
eral class of gas flows: 

Consider gas flows in which (A) the pressure differences arise 
through inertial effects and (B) the Mach Number is reduced by 
reducing the flow velocities. From condition (A), the variations 
with velocity (u, say) are 


p(Ou;/Ox;) ~ u (since p is fixed) | 

Dp Op (4) 
~~ Ub; ~ u (since Ap ~ u?*) \ 

Dt Ox; 

It follows that, for such low-speed gas flows, the energy equa- 
tion is better approximated by 

pc)y(DT/Dt) = (0/dxx) {k(OT/dx,)} (5) 
than by Eq. (2). In such cases, p(Ou;/Ox;) cannot be directly 
neglected in Eq. (1). 

Pressure differences arise through inertial effects, for example, 
in all flows having a body immersed in an extended fluid flow with 
p = constant at infinity in all directions. The fluctuation field in 
homogeneous turbulent gas flow also falls in this ‘‘inertial’’ class, 
and the energy equation in reference 2 therefore applies to low- 
speed gas flow. For liquids, c, replaces Cp. 

For gas flows with the pressure differences imposed from out- 
side (e.g., Poiseuille motion), we have reached no general con- 


«clusion. 
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Comments on the ‘‘Effect of Surface Cooling 
on Laminar Boundary-Layer Stability’’ 


Lester Lees 

Aeronautical Engineering Laboratory, Princeton University, Prince- 
ton, 

September 21, 1951 


I REPLY TO MARTIN BLOOM’s NOTE on the effect of surface 
cooling on laminar boundary-layer stability (Readers’ Forum, 


AL SCIENCES—DECEMBER, 1951 


JOURNAL OF THE AERONAUTICAL SCIENCES, September, 1951), 
his conjecture about a numerical error in the calculations of min- 
imum critical Reynolds Number at high supersonic Mach Num- 
bers in N.A.C.A. T.R. No. 876 is correct; the temperature gra 
dient was calculated incorrectly in these cases. Fortunately, the 
calculation of the effect of cooling at J = 0.70 and the caley- 
lations for the insulated wall up to M = 1.3 are not affected. 

van Driest’s work, published in the Readers’ Forum, JourRNAL 
OF THE AERONAUTICAL SCIENCES, October, 1951, shows trends 
similar to Bloom's, but the numerical values do not agree, proba- 
bly because van Driest utilized the more exact Sutherland viscos- 
ity-temperature relation. The physical basis for the instability 
at high supersonic speeds is the large positive density gradient in 
the outermost portion of the boundary layer. Because of this 
gradient, the important quantity (d/dy)[p(du/dy)| is also posi 
tive up to large values of wave velocity. 

It would not be fair to draw any conclusions as yet as to transi 
tion at high supersonic speeds, at least until calculations or theo- 
retical estimates are available of the rate of amplification of the 
unstable disturbances, and the role of the nonlinear effects in 


transition 


On the Flow of Viscous Incompressible Fluids 
in Prismatic Ducts 


Martin Lessen 

Professor, Department of Aeronautical Engineering, The Pennsylvania 
State College 

September 21, 1951 


_—o FLOW OF A VISCOUS INCOMPRESSIBLE FLUID through a 
prismatic duct can be expressed as follows: 


ou (> a) 
p = P(t)h+u + (1) 
ol oy? 02? 
where 
P(t) = —vup/ox (function of time) 
x = dimension in direction of flow 
¥, 3 = dimensions perpendicular to direction of flow 
t = time dimension 
p = pressure 
p = density 
m = absolute viscosity 
u = velocity of flow 
First consider the case where P(t) = 1 and uv = 0 when? = 7 
Of course, « = 0 also at the boundary at all times. For/ > 71, 


the solution will be the same as the solution for a unit step func 


tion in pressure gradient with the flow starting from rest at ¢ = 
r. Fort >+r,u willbe of the form 


u= F(y,2,t — 7) (2 


The general solution of Eq. (1) therefore is of the form 
J t , . 

u(y, c, t) = PoF(y, zs, t) + 5 P'(r)F(y,s,t —7r)dr (3 
where P,) is the pressure gradient at the beginning of events, 
(¢ = O)and uv = 0 when? < 0. 

Consider the flow in a rectangular duct of dimensions 2a by 
2b such that 
u = O0aty = ta 
u= Oats = +b 


It follows that 








an 


all 
th 


ch 
of 

co 
wi 
ed 


tu! 


no 


an 


Ser 


Sep 





1951), 


f min- 


Num- 
e gra 

V, the 
-aleu- 


"RNAL 
rends 
roba- 
iscos- 
bility 
“nt in 
F this 


posi 


ransi- 
theo- 
f the 
ts in 


ds 


yania 


gh a 


nts, 


. by 





READERS’ 


Fly, 2,¢ — 7) = 
— : OE 
64a2b2 >» = (—1)[On+n) 2 1 
ter’ ‘ - 4 homed. mn( mb? + na?) 
‘ m=1,5,0 n=1,9,0 
mnry NTs > 9 > 2 »\ 
cos =i cos ah [ — et p) [(m2x2/ 4a?) + (n2x2/4b2) (tf —7 ] (4) 
2a 20 


Eq. (4) inserted into Eq. (5) will give the general solution of 
u for a time variable P for a rectangular prismatic duct 


A Note on Wind-Tunnel Wall Effect Corrections 
for Rectangular Airfoils 


K. G. Merriam 
Professor, Aeromechanics, Worcester Polytechnic Institute 
August 16, 1951 


a A WIND-TUNNEL TEST of an airfoil, one of the re- 
quirements is to compute the rate of change of lift coefficient 
with respect to angle of attack, corrected to infinite aspect ratio 
and free air conditions. 

In such a computation, a useful reference fact is that, for an 


airfoil of infinite aspect ratio at angles of attack below burble, - 


the value of drag coefficient is essentially constant. Therefore, 
it has essentially no rate of change with respect to anything, in- 
cluding the square of lift coefficient. 

To show how this fact may be used, call A the value of the rate 
of change of drag coefficient with respect to the square of lift 
coefficient, as found from the slope of the plot of the appropriate 
wind-tunnel data. Define 6 and 7 as on page 147 of the 1948 
edition of Elements of Aerofoil and Airscrew Theory by Glauert 
Also, define € as a correction term that exists because of wind 
tunnel wall effect. 

In the light of the reference fact mentioned above, one may 


now write 
A — }[(1 + 6)/r A.R.] +e} =0 


and solve for « as the single unknown. 

Now define B as the rate of change of angle of attack, in de- 
grees, with respect to lift coefficient, as found from the slope of a 
second plot of appropriate wind-tunnel data. 

One can now express the required rate of change of lift co- 
efficient with respect to angle of attack, in degrees, corrected to 
infinite aspect ratio and free air conditions, as the reciprocal of 
the quantity: 


B — 57.3} [(1 + 7r)/w A.R.] + ef 


In the absence of more reliable values of 6 and 7, the theo 
retical values given in Fig. 85, page 147, of the Glauert reference 
mentioned above should lead to only minor error. 


Comments on ‘The’ Laminar-Turbulent 
Transition in a Boundary Layer’’ 


G. S. Hislop 
Senior Assistant to the Controller of Research and Special Develop- 
ments, British European Airways 


September 20, 1951 


IT Mr. EMMons’ PAPER, ‘‘The Laminar-Turbulent Transition in 


a Boundary Layer—Part 1,”’ given in the July, 1951, JOURNAL, 


FORUM 845 


he refers (page 491) to the “‘spot’’ of turbulence growing at each 
point normal to its surface by consuming the laminar region 
TFiom this it might be inferred that at a given downstream sta- 
tion on a flat plate, such as was employed in his experiments, the 
proportion of turbulent to laminar flow over a period would de- 
crease with distance from the surface of the plate. 

Evidence that this is so was found during the course of experi 
ments at Cambridge University in 1939 and described in some 
unpublished work.! 

For example, using a 2!'/2/10,000 diameter hot wire 5*/,-in. 
downstream from an 0.020-in. high ridge affixed transversely 
across the surface of a glass plate, at a wind speed of 110 ft. per 
sec., and a Reynolds Number of transition of 2.4 X 10°, the flow 
at 0.007 in. from the surface was turbulent for approximately 
60 per cent of the time and at 0.050 in. for roughly 40 per cent of 
the time 

Mechanical trouble with the oscillograph camera prevented 
me from obtaining a satisfactory record at the time. However, 
there was no doubt of the diminishing proportion of turbulence 
as the hot wire was withdrawn from the surface of the plate. 
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Uniformly Loaded Semi-Infinite Wedge-Shaped 
Plates 


A. J. A. Morgan 
Aeronautical Engineering Research Consultants, Pasadena, Calif 
September 14, 1951 


hes DEFLECTIONS OF A UNIFORMLY LOADED semi-infinite wedge- 
shaped plate of constant thickness (Fig. 1) are explicitly de- 
termined for the case where one edge is free and the other is 
clamped. 

Using the notation of Timoshenko,' the deflections w(x, y) of a 
constant thickness plate are governed, under the assumption of 
small slopes, by the partial differential equation 

4ay Ow Ow P(x, y) 
Vw = +2 + = : (1) 
ox! Ox*dy? oy! D 
where P(x, y) is the loading (force/unit area) and D is the plate 
stiffness factor. 
The form of Eq. (1) is unchanged under the transformations 


x¥=a"x, Vra"y, w = arw (2)* 
where a(-+£0) is a numerical parameter and m, n, p are real con- 
stants, if 
p/m = 4 (3) 
and 

P(X, 3) = P(x, y) (4) 


The condition (4) is satisfied, for example, when the load is con 


stant. Two quantities that remain unchanged (absolutely in 
variant ) under the transformations (2) are 


ai"™, w/xPp/m (5) 
Hence, by (3) and (5), the substitution of 
w= x'f(n), n = y/x (6) 
into Eq. (1) will reduce it to the ordinary (similarity) differential 
equation? 
* Mathematically, this is a continuous one-parameter group of trans 


formations. 
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2 d‘f e d*f The boundary conditions at the clamped edge are 
wo > iP = > Se + $) —. + 
dn dn w=0 (zero deflection ) (9) 
(Sn? +1 d*f 4 df 4 24f P (7) and 
4(3n* ) — 24 <4f = . 
ai dn? dn J ps “ 
w ow ow 
The general solution of Eq. (7) is ae lox tan B — >y | °° 8 =0 (zero slope) (10) 
Nah a Tee te ee (1° _ L Cn? + Cn* — 39%) + LP 7? where n is the outward pointing normal to the edge of the plate 
3 & D rhe boundary conditions at the free edge’ can, for this case, be 
(8) expressed as 


‘ Ww " Ow , Ow . 
ww + (1 — v) - cos? B + = au f= sin 28 | = O 11) 
Ox? oy* Oxoy 
and 
ee o .. 3) O*w 1 {/o*w O*w 
Vw sin B + V+w cos 8 (1 v) : cos 28 + = 0 12 
Ox oy Os LOxoy 2 \oy? Ox? 
where 
V2w = (0%w/Ox?) + (07w/dy?) 
s = tangential direction at free edge (positive clockwise) 


0/05 = —cos B(0/dx) — sin B(d/dy) 


The conditions (11) and (12) correspond, respectively, to the requirements that the bending moment and shear force must vanish at 
the free edge. 

Consider the clamped and free edges of the plate to be placed at n = —k = tan(—8)andy = k = tan 8, respectively. The condi- 
tions (9) to (12) can, by Eqs. (6), then be expressed as 


f(-—k) =0. (13) 
(df/dn), = —-~ = 0 (14) 
ee df : d*f 
12(1 + vk?)f(k) — 6k [vk? + (2 — v)] . + [vk* + 2(2 — v)k? + v] we = 0 (15 
dn] y =k dn? ] =k 


and 


l-— yp df d*f d*f 
24k (1 + f(k) + 6[8k? + (2 — v»)][ - — 6(1 + k?) + (1 + k?2)4 — = 0 (16 
1 + k* dn/n=k dy? ] 4 =k dn *],=k 


The expressions (13) to (16) are entirely in terms of f() and its derivatives; hence, the problem can be completely solved in terms of 
the solution (8). 

As a simple example, consider the clamped edge of the plate to be placed at the x-axis and the angle 8 to be such that 0 < 6 < 7/2 
On altering (13) and (14) accordingly, (8) becomes: 


S(n) = Cyn? + Ci(n*t — 3?) + (1/8) (P/D)n? (17) 


The substitution of Eq. (17) into Eqs. (15) and (16) gives two simultaneous equations for C; and Cy—namely, 


LP 
6vk(k? + 1)C3 + 2[k4(6 — v) + 2k1 + v) — odCGy = dD [2k2(1 — 2p) vki vy] (18 
and 
2 peer : " : kP 
8k® + 7RUZ — v) + 3RA1 — v) +: 1)C; + 2k[k5 — 38vy) + (1 + »)IG = 1D [4k* + k(5 vy) + 2k(1 —»)+(1—» (19 


Hence, in determinant form, 


Er 4 2k2(1 — 2v) — vk! v k4(6 vy) + 2k*(1 + v) v on 
= vo _ : oe <U) 
2D A —|4ke +4 R(5 - 9) + 2k(1 —yv)+k1 —v)] kRY5 — 38v) + k1 +7 
and 
C Te 6vk(k? + 1) 2k(1 2v) — vk‘ v 21 
; 4D A! 8k® + 7k4(2 vy) + dR7>(1 vy) + 1 [4k5 + R95 v) + 2k%(1 v) + Rk(1 — pv) ia 
where 
oo 6vk(k? + 1) R46 v) + 2k%(1 + v) v 


8k® + 7RU2 — v) + 3R%~1 v) +1 R35 3v) + Rk(1 + v) 





=<=— —— 


Tht 
the 


cor! 


suri 


ana 


whi 


Th 
str 
we 


we 


thi: 


and 





(9) 


) (10) 


plate 


ase, be 


(11) 


(12 


ish at 


ondi- 


(16 


1S of 


(17) 


19 


20 ) 





READERS’ 











The semi-infinite wedge-shaped plate 


Fic. 1. 


Thus it is seen that the variation of the deflection pattern with 
the angle 8 = tan~! k obeys a complex law. 


Consider again the case shown in Fig. 1. The stresses at the 


corner (x = 0) shall be determined. The stresses at the upper 
surface of the plate are given by* 
Et o*w o*w 
" ; tii Pe zs Me ak . 
(1 — »*) \ Ox? oy" 
Et Ow dow 
eC : — 7. ee 
(1 v*) \Oy- Ox? 
and 
Tsy ™ - [Et (1 , 2 )] (0?w Oxoy ) 


which in terms of f(7) and its derivatives can be expressed-as: 


Etx? ; df d*f 
12f(n) — 6n — + (9? + ») =, (22) 
(1 — py?) dn dn? 


1: = 
Etx? ; . df _. oF a 
vy = 12vf(m) — 6un t (1 + °°) —. (23) 

(1 — pv?) dyn dn? 
Etx? 3 df =) 04 
Tzy ™ . U] . (24) 

(1 + pv) dn dn? 

The above expressions show that, on any line » = constant, the 


stresses will depend only on x. Hence, as x — 0 (the corner), 
we have oz, oy, Try > O at the same rate that x* — O for any 
wedge angle 28 less than 180°.* 


* Plates with a wedge angle greater than 180° cannot be considered by 
this method since there is no mathematical distinction between the right 
and left-hand sides of the (x, y)-plane 
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Although the expressions given above constitute a mathemat- 
ically exact solution, under the physical limitations of the theory, 
to the problem posed, it is evident that the solution has physical 
significance only for those values of x for which the stresses do not 
exceed the proportional limits of the material [see Eqs. (22) 
(24)]. Let 6 denote the least such x. 

Thus, Eq. (8) can be considered as the exact solution for finite 
triangular plates with a particular stress distribution, Eqs. (22)- 
(24), at an arbitrary crosswise edge, 

(25) 


y = f(x), 0O<x<b 


placed as shown in Fig. 1. In particular, (25) can be the straight 
line x = bor the arc of a circle 

x? = §?/(1 + 7?) (26) 
in which case it is called a sector plate. 

The stress distributions, Eqs. (22)-(24), can then be looked 
upon as boundary conditions at the crosswise edge, Eq. (25), of a 
triangular plate, and they are such that no stress singularity 
(or, Gy, Try — ©) occurs at the corner for any wedge angle 0 < 
28 < (Fig. 1). 

The deflection 
wedge-shaped plates satisfying different conditions at the radial 
edges (clamped-clamped, clamped-simply supported, simply sup 
ported-simply supported, etc.) can also be determined by the 
method outlined above, since all of the boundary conditions can be 
expressed in terms of f(n) and its derivatives 
tained in such a manner can also be considered to give, under the 
provisions outlined in the preceding paragraphs, the deflections of 


patterns for uniformly loaded semi-infinite 


The solutions ob 


finite triangular plates. 

It should be noted that the method is also applicable when the 
load is not uniform but is a function of y alone [that is, P(x, y) is 
invariant under the group (2)]. Especially tractable will be the 
case where P(x, y) is a polynomial in », since the particular solu- 
tions of Eq. (7) can then be obtained with ease. 

The author is indebted to Prof. A. D. Michal of the Depart- 
ment of Mathematics, California Institute of Technology, for the 
suggestions extended during the preparation of this note. 
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